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We present a new evolution system for Einstein's field equations which is based on tetrad fields and 
conformally compactified hyperboloidal spatial hypersurfaces which reach future null infinity. The 
boost freedom in the choice of the tetrad is fixed by requiring that its timelike leg be orthogonal to 
the foliation, which consists of constant mean curvature slices. The rotational freedom in the tetrad 
is fixed by the 3D Nester gauge. With these conditions, the field equations reduce naturally to a 
first-order constrained symmetric hyperbolic evolution system which is coupled to elliptic equations 
' for the gauge variables. The conformally rescaled equations are given explicitly, and their regularity 

at future null infinity is discussed. Our formulation is potentially useful for high accuracy numerical 
modeling of gravitational radiation emitted by inspiraling and merging black hole binaries and other 

■ highly relativistic isolated systems. 
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The calculation of gravitational waveforms generated by the inspiral and coalescence of binary black holes is now a 
routine matter for numerical relativists (see Refs. [H for recent reviews). One of the most pressing goals for those 
working in the field today is to compute these waveforms as accurately as possible, without losing computational 
efficiency. Optimally, this goal is accomplished by including future null infinity (^ + ) in the computational grid 
because then the Bondi news function [3j-l5|, which contains all the gravitational wave information, can be read off 
without the need for extrapolation. One way to include J? + in the numerical grid is Cauchy characteristic matching 
(see Ref. for a review), which joins a standard 3 + 1 Cauchy code in an interior region with a characteristic code in 

■ an outer region. However, arranging stable data transfer in both directions at the interface between the two regions 
' is difficult in general; hence, a modified approach called Cauchy characteristic extraction (CCE) , m which data 

is extracted from the interior of a standalone Cauchy development and extrapolated to using a characteristic 
code, is now often used to extract binary black hole inspiral and merger waveforms at . While CCE is certainly a 
significant improvement over wave extraction at the finite boundary of a conventional Cauchy code, potential sources 
of error remain. Initial conditions for the characteristic code cannot easily be made compatible with the Cauchy 
time development. Also, either errors in boundary conditions for the Cauchy code can affect the CCE or there is 
insufficient time for all the junk radiation present in the initial data of the Cauchy code to escape. Furthermore, 
the extraction imposes a gauge on the characteristic code in which the expansion of generally does not vanish, 
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I. INTRODUCTION 



X 

^ . making it somewhat of a challenge to obtain the Bondi news accurately [8[ . 

An approach which, we argue, has greater potential for high accuracy waveform calculations is to conformally 
compactify the spatial geometry on hyperboloidal spacelike hypersurfaces, as originally advocated by Friedrich but 
with elliptic equations determining the gauge degrees of freedom. Hyperboloidal hypersurfaces behave like conventional 
time slices near the compact object sources, but smoothly transition to asymptotically approach null slices in the 
physical spacetime at . For a wide class of problems of physical interest, the compactified conformal geometry is 
regular, and the entire space section can be represented on a finite resolution and relatively small coordinate grid, 
enabling (we expect) efficient simulations and accurate waveform calculations. We can choose boundary conditions 
on the elliptic gauge equations at ^ + which make the expansion of vanish, and thereby simplify reading off the 
asymptotic gravitational wave amplitudes. No outer boundary conditions are required on the dynamical variables, 
since ^ + is an ingoing null hypersurface in the conformally compactified geometry. 

Friedrich Q derived a first-order symmetric fully hyperbolic system evolved on compactified hyperboloidal hyper- 
surfaces, with tetrad components of the Weyl tensor as the fundamental dynamical variables and with the Bianchi 
identities providing the evolution equations for these variables. The system is manifestly regular at J? + ', thanks to 
the conformal invariance of the Weyl tensor. We refer the reader to the review article [Io| by Frauendiener for further 
details, including references to numerical work based on this formulation. The system is very useful for obtaining 
analytic results on existence and uniqueness of radiative spacetimes (e.g. fill]). 

Metric-based proposals for solving the Einstein equations on conformally compactified hyperboloidal hypersurfaces 



2 



include the completely hyperbolic generalized harmonic scheme of Zcnginoglu [12j and the mixed hyperbolic-elliptic 
formulation of Moncrief and Rinne Both of these are second order in spatial derivatives. Rinne [l4[ has imple- 
mented a simplified version of the Moncrief-Rinnc formalism in an axisymmetric finite difference code, and demon- 
strated long-term stable dynamical evolution of the Einstein equations. Zcnginoglu and Tiglio fl5j and Zenginoglu 
and Kidder [l6j have successfully evolved wave equations for test fields on compactified hyperboloidal hypersurfaces 
on fixed black hole backgrounds. Conformally flat hyperboloidal Bowen-York binary black hole initial data have been 



constructed in 17 1 



The tetrad framework we present in this paper uses the dyadic notation of the Einstcin-Bianchi system of Estabrook 
and Wahlquist [l8j] , modified to allow for conformal rescaling. Instead of basing the dynamical system on the Weyl 
tensor and the Bianchi identities, our fundamental variables are the 24 tetrad connection coefficients (i.e., Ricci 
rotation coefficients) evolved using the Einstein equations. In the context of certain dynamic al g auge conditions, this 
system can be put in a first-order symmetric hy per bolic form, the "WEBB" equations [ill [2(| ■ The mathematical 
structure of these equations is analyzed in Refs. j2lll22l|. and ID numerical tests are presented in Ref. (23|. A more 
general discussion of symmetric hyperbolic systems for the Einstein equations and their suitability for numerical 
relativity is given by Friedrich in |24|. Purely hyperbolic numerical evolution in general relativity, particularly in a 
first-order scheme, needs to accurately preserve a fairly large number of constraints. 

Our approach here, like that of Moncrief and Rinne (l3j . is a mixed hyperbolic-elliptic system based on a special 
class of hyperboloidal hypersurfaces called constant mean curvature (CMC) hypersurfaces, for which the trace of the 
extrinsic curvature is uniform in both space and time. Some of the elliptic equations are singular at J* + , where the 
conformal factor vanishes, which forces the solutions to have a highly constrained asymptotic behavior. This is a 
feature, not a bug, since it enables a detailed analysis of the asymptotic regularity conditions which must be imposed 
in order to obtain non-singular evolution of the conformal geometry. We differ from Moncrief and Rinne in that 
our evolution equations arc first order in space as well as time, and include connection coefficients of the conformal 
spatial geometry in our fundamental set of dynamical variables. Also, instead of the coordinate metric and coordinate 
components of tensors, our approach is based on orthonormal tetrad basis vectors. The tetrad components of all 
tensors are coordinate scalars, and the metric for raising and lowering tetrad indices is trivial. Consequently, the 
resulting equations are simpler in form and in interpretation than those of Moncrief and Rinne. 

The extra gauge degrees of freedom in our tetrad formalism corresponding to local Lorentz boosts and rotations 
of the tetrad frames are fixed in a natural way which simplifies the equations and reduces the number of variables 
compared to a general tetrad formalism. The time-like vector of the tetrad is chosen to be the unit normal to the CMC 
hypersurfaces, which constrains the triad of spatial vectors to be tangent to the CMC hypersurfaces. The rotational 
degrees of freedom in the conformal spatial triad vectors are fixed by imposing the 3D Nester gauge [25[ . In a certain 
global sense, the Nester gauge makes the triad vectors at different points in the hypersurface as nearly related by 
parallel transport as possible. The Nester gauge leads naturally (2(| to a choice of conformal factor, which satisfies 
an elliptic equation derived from the Hamiltonian constraint and is enforced by elliptic equations for potentials that 
determine the conformal angular velocity of the triad frames and the trace of the conformal extrinsic curvature. Of 
the 24 tetrad connection coefficients, only 10 are dynamical and form a Maxwell-like symmetric hyperbolic system. 
The remaining evolution equations are the nine linear advection equations for the coordinate components of the three 
spatial triad vectors. 

As in Ref. fl3| . we put considerable effort into analyzing the properties of the equations and their solutions in the 
neighborhood of . The evolution equations for the dynamical variables contain terms which have a conformal 
factor in the denominator and therefore are potentially singular at ^ + . Certain well-known regularity conditions 
must be imposed to make these terms finite, such as the "zero-shear" condition on the null generators of <f + (see 
Ref. [13]). Once imposed in the initial conditions, they are preserved by the evolution equations. An additional and 
stronger regularity condition is also usually assumed, that the conformally invariant Weyl tensor vanish at J? + . We 
refer to this as the "Penrose reg ularity condition", and it is also automatically preserved by the evolution equations. 
The review of Frauendiener [lj| puts the various regularity and smoothness assumptions at J zr+ in context. 

The question of how much smoothness in the conformal geometry at .y + is, in some sense, physically generic has 
been the matter of considerable debate in the literature. The original analyses of Bondi et al. H and Sachs [H 
assumed analyticity, but the more systematic conformal analysis of Penrose [28j reduced this to C 3 smoothness. 
Later, Chrusciel et al. [29| showed that "polyhomogeneous" terms in the expansion of the angular part of the metric 
in Bondi-Sachs coordinates, terms of the form x n (logx) m , where x = 1/r, are consistent with the Einstein equations 
and the ability to define a physically suitable Bondi energy. The coefficient of the "leading log" term is a constant 
during evolution, so if these polyhomogeneous terms are ever present they were always present back to past null 
infinity, where they are associated with incoming radiation. 

We conclude from this that in the context of typical astrophysical problems, where one can assume there is no or 
negligible incoming radiation at past null infinity, it is perfectly appropriate to exclude polyhomogeneous terms of the 
type contemplated in Ref. [29] from ever being present. Similarly, violations of the regularity conditions described 
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above, including the Penrose regularity condition, can only be present if they were present at past null infinity, and 
cannot be generated by any dynamics of the physical system. However, as was shown in Ref. [27j and as we derive 
in Appendix [A] polyhomogeneous terms in both the conformal and physical extrinsic curvature and the conformal 
factor are generically present when the spacetime is foliated on CMC hypersurfaces. In fact, we show that they must 
be present whenever outgoing radiation is present at , and are not present otherwise. These terms arise from the 
special character of the momentum constraint equation and the elliptic equation for the conformal factor on a CMC 
slice as it becomes asymptotically null at J^ + . In a more flexible hyperboloidal slicing, where the trace of the extrinsic 
curvature only approaches a constant at , the polyhomogcous terms need not be present. In fact, they are clearly 
not present when the spacetime is evolved using the regular conformal field equations of Friedrich 0] , which allows 
the foliation to evolve as part of the hyperbolic system. 

There are some significant advantages to our mixed hyperbolic-elliptic system. The boundary conditions on the 
gauge variables determined by the elliptic equations can be adjusted to simplify and control the geometric properties 
of and coordinates on J^ + . For instance, the boundary condition on the elliptic equation for the conformal lapse, 
together with the boundary condition on the trace of the conformal extrinsic curvature, can keep the intrinsic conformal 
geometry of the intersection of the CMC hypersurface with J? + that of a two-sphere of constant area and the time 
coordinate on the CMC hypersurfaces equal to the retarded Minkowski time at J^ + . The shift vector controls the 
evolution of the spatial coordinates. We propose a few alternative vector elliptic equations for the shift, all with a 
boundary condition chosen to keep at a fixed coordinate radius. 

The final conformally rescaled system is a symmetric hyperbolic evolution scheme which is coupled to an elliptic 
system for the conformal factor and the gauge variables. Hyperbolic-elliptic formulations have been discussed in 
the literature, both from mathematical and numerical viewpoints. A hyperbolic-elliptic formulation of Einstein's 
equations in spatial harmonic coordinates on CMC hypersurfaces, very similar to the Moncrief-Rinne formulation 
but without the conformal compactification, has been proven well-posed by Andersson and Moncrief (30| . Note that 
the Andcrsson-Moncrief formulation is not constrained - elliptic equations are solved only for the lapse and shift. 
The Meudon group proposed a fully constrained, conformally rescaled system [3l| for which both numerical tests and 
mathematical analyses have been performed (3ll - l33| . Finally, elliptic equations have been periodically enforced in a 
constraint projection scheme 3J] for a first-order reduction of a scalar field equation on a black hole background, 
and shown not to be prohibitively expensive. Therefore, we hope that a numerical implementation of the prop osal 
presented in this paper with efficient elliptic solvers such as those in the Spectral Einstein Code SpEC [3^, |36[ will 
be possible with reasonable computational effort, especially since we would need relatively small grid sizes due to 
conformal compactification and the relatively small number of wave oscillations out to J? + on a CMC hypersurface. 

Our paper is organized as follows. In Sec. [Hi we specify our choice of variables and notation and, based on the 
hypersurface-orthogonal gauge, we present the 3 + 1 decomposition of the triad and connection coefficients as originally 
developed in the dyadic formalism of Ref. [l8j]. In Sec. IIII1 we discuss the evolution and constraint equations of the 
WEBB tetrad formulation when specialized to a hypersurface-orthogonal gauge, and we decompose both the extrinsic 
curvature and the dyadic form of the spatial connection coefficients into symmetric trace-free, antisymmetric, and 
trace parts. The antisymmetric and trace parts are fixed by our gauge conditions which are discussed in Sec. IIVI 
While the antisymmetric part of the extrinsic curvature is identically zero, its trace is fixed to a constant by the CMC 
slicing condition, and the trace and antisymmetric parts of the dyadic form of the spatial connection coefficients are 
determined by the 3D Nester gauge condition (the latter up to a gradient). The remaining symmetric trace- free parts 
are governed by a nice Maxwell-like symmetric hyperbolic evolution system which is coupled to advection equations 
for the triad vectors and elliptic equations for the gauge variables. The conformal decomposition is performed in 
Sec. |Vj The Nester gauge is covariant with respect to conformal transformations, and this property motivates a 
natural choice for the conformal factor. With this choice, the antisymmetric part of the dyadic representation of the 
conformal connection coefficients is zero, which further simplifies the equations. Also in Sec. |V1 we derive manifestly 
regular elliptic equations guaranteeing that these gauge conditions are preserved throughout evolution, and we discuss 
possible elliptic equations for the shift. In Sec. IVI[ we analyze necessary conditions to impose on the fields at J? + in 
order to have a regular time evolution. Furthermore, we evaluate explicitly the apparently singular terms at ^ + in the 
right-hand side of the evolution equations. Our evolution system, together with its constraints and elliptic equations 
for the gauge variables including boundary conditions, is summarized in Sec. I VIII and issues such as comparisons 
with other approaches and extraction of gravitational waveforms are discussed in Sec. I Villi More technical aspects 
relevant to our formulation are given in Appendices E1[B] and [C] where we discuss, respectively, asymptotic expansions 
of the fields near , the constraint propagation system, and the relation of the 3D Nester gauge condition to the 
3D Dirac equation. 

Throughout this paper, we use the sign convention (— , +, +, +) for the metric. We denote by the 2D intersection 
of the .y + null hypersurface with a hyp erboloidal CMC spatial hypersurface, and by the = sign an equality that holds 
only at ,y + . We adopt the Wald [331 s i§ n convention for the extrinsic curvature of a spatial hypersurface, so it is 
positive for expansion toward the future. 
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II. 3 + 1 DECOMPOSITION 



The tetrads in the spinor formulation of Newman and Penrose 38[ are aligned along null congruences and are 
therefore conducive to understanding radiative properties of a spacetime. In contrast, in the tetrad formulation 
presented by Estabrook and Wahlquist [Tsj . the tetrad is oriented with respect to a preferred timelike congruence. 
This can be useful if a physically preferred timelike direction, such as a fluid 4-vclocity, exists. However, this is not 
the case in vacuum. Our choice is to orient the timelike leg of the tetrad orthogonal to the CMC hypersurfaces of 
constant coordinate time S t . It follows then that the spacclikc legs of the tetrad are tangent to S t . This is called a 
hypersurf ace- orthogonal gauge condition, and we assume this gauge condition throughout the paper. 

Our notation follows that of Ref. [HI]. Spacetime indices (0 — 3) are Greek and spatial indices (1 — 3) are Latin. 
The beginning of the alphabet is used for tetrad indices, and the middle of the alphabet is used for coordinate indices. 
Spatial tetrad indices will not be raised and lowered, since the spatial metric in the tetrad basis is trivial. In Ref. [T^j, 
there was a distinction between the tetrad vectors with components \£ and the spatial triad vectors B a defined as 
projections of the spatial tetrad vectors tangent to S t , with only spatial components . In a hypersurface-orthogonal 
gauge, the spatial tetrad vectors have zero time components and are identical to the B a . The directional derivatives 
along the three spatial legs of the tetrad are then 1 



The Lie derivative along the timelike leg of the tetrad is written in terms of the lapse function (a) and the shift vector 

(/3 fe ): 

Do = -(i-£ ). (2) 



a \dt 

Let us stress that we regard variables with pure tetrad indices (such as K a b, N a b, a& and Ub below) as scalar fields on 
St. In this sense, the Lie derivative operator £p with respect to the shift acts as the directional derivative on such 
variables. In contrast to this, we have olDqB* = (dt — ftd^B* + D a fi % which differs from the definition in [l^] , where 
aD B: = (dt~^d 3 )B a \ 

The coordinate basis spacetime metric g^ v and the spatial metric hij can be calculated from the tetrad and triad 
vectors as 

g r = rfl i \»\g i h«=B*B>, (3) 

where rj a ^ is the Minkowski metric. The connection coefficients with respect to the Levi-Civita connection V of g^ 
in the tetrad basis are calculated from the commutators, 

T a fj 7 = X a ■ V 7 A^j = — r ( g Q7 = - (A^g • [A Q , A 7 ] + A 7 • [A Q , \p] — \ a ■ [Xp, A 7 ]) . (4) 

In the dyadic formalism (l8| . the 24 independent T a p 7 are decomposed into two 3x3 matrices 

Kab = TfeOa, N a b = - Ebcd ^cda, (5) 

and two vectors 

a-b = Tsoo, ^b = — - Ebcd L c do- (6) 

In the hypersurface-orthogonal gauge, the K a b are the triad components of the symmetric extrinsic curvature tensor 
of Sf. The sign of K a b is opposite to that of Misner, Thorne, and Wheeler (39|, in that positive values imply expansion 
of the normals, not contraction. The N a b describe the induced connection coefficients of S t . It is often convenient to 
write the anti-symmetric part of N a b as the vector 

nb = ^ e h cd N cd . (7) 

The vectors ab and Wf, are, respectively, the acceleration and the angular velocity relative to Fermi- Walker transport 
of the tetrad frame. The acceleration vector field is given by the gradient of the logarithm of the lapse, 

a b = D b (loga). (8) 



1 The vector A a of Ref. [Til , which is the velocity of the tetrad congruence relative to a hypersurface rest frame, is identically zero. 
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III. BASIC EQUATIONS 

Here we write the basic Einstein evolution and constraint equations of the WEBB tetrad formulation (these are 
derived in Ref. [l9| from the Einstein equations and the Riemann constraints, see also Appendix [BJ, specialized to 
a hypcrsurfacc-orthogonal gauge. Note that K ab is explicitly symmetric. The resulting evolution equations for K ab 
and N ab , respectively, arc 



D K ab - e a cdD c N db = (D a + a a )a b - e bcd N ac a d + 2 e cd ^ a K b y uj d + NN ab 

1 



D N ab + e acd D c K 



ill, 



+ ^e adf e bce (K dc K fe - N dc N fe ) - K ac K cb - N ca N cb + 8nG 

= -(D a + a a )u> b + e bcd (K ac a d + N ac u d ) + e acd N cb uj d 
- NK ab + 2N c [ a K b ] c + e ad fe bce N dc Kf e — 8irGe abc j c , 



Oab - -^5 ab (a cc + p) 



(9) 



(10) 



where N is the trace of N ab , K is the trace of K ab , and a b = D b (\oga). Here, we have also defined the energy density 
P = T 00i the energy flux j b = —rob, and the stress tensor a ab = r abl where t^ v denotes the stress-energy tensor. 
The Hamiltonian and momentum constraint equations are, respectively, 



2D a n a = N ab N ab + - (K ab K ab - A^ ab A^ ba - K 2 - iV 2 ) + 8irGp, 



D h K 



b^ab 



D a K 



Saba K bd N dc + 2K ac n c - 8irGj a . 



(11) 
(12) 



There is also a constraint for N ab analogous to the momentum constraint equation, which is obtained from the 
Riemann identities: 



D b N ab - D a N 



^abc N bd N dc . 



(13) 



Evolution and constraint equations for the spatial triad vectors B a are derived from the commutators of the 
orthonormal basis vectors and are 



D B a k 

£abc D b B 



—e abc oj b B c — K ac B c , 
N ba B b k -NB a k . 



(14) 
(15) 



In the following, it will be useful to separate out the symmetric, trace-free parts of K ab and N ab , denoted by K ab 
and N ab , so that 

Kab = K ab + ^ K, N ab = N ab + £ab c n c + ^j-N. 



From Eqs. (j9|) and (fTQ|). we obtain the following evolution equations for K ab and N ab : 



2N ~ 

D K ab - D c N d(a e b)cd = <{ (D (a + a {a )a b) - {D {a - a {a )n b) - K K ab + —N a b - 2N ac N cb 



z -cd(a 



N b)c (2n d - a d ) + 2K b)c ui d 



SttG, 



TF 

Cab J" 



(16) 



D N ab + D c K d ( a £b) 



cd 



AN - K - - 



~'cd(a 



K b)c a d + 2N b)c Lo d \ (17) 



TF 



where {...} TF denotes the trace-free part of the expression inside the parenthesis. The evolution equation for the 
trace of the extrinsic curvature, from which an elliptic equation for the lapse will be derived in Sec. IIV A[ is 



D a K = (D b + a b - 2n b )a b - K ab K ab - -K 2 - AitG((j aa + p). 



(18) 



We have used the Hamiltonian constraint, Eq. (jllj) . to eliminate the scalar 3-curvature. The evolution equations 
for the trace and the antisymmetric part of N ab , from which we will derive elliptic equations for u a and the time 
derivative of the conformal factor in Sec. IV Al are 



D N = -(D a + a a )uJa- N ab K ab - -KN, 

2K 



(19) 



2D a n a = -s abc (D b + a b - 2n b ) u c — (a a + n a ) - e abc K bd N dc + K ab {a b + n b ) - 8-KGj a - (20) 
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The second equation has been simplified by applying the momentum constraint equation (|12[) . 

The Hamiltonian and momentum constraint equations and (|12p , expressed in terms of the symmetric trace- free 
and antisymmetric parts of K a b and N ab , are 

\ / - - - - \ 3 (K 2 + N 2 ) 
D a n a = - [K ab K ab + N ab N ab j + - n a n a - ± '- + AnGp, (21) 

D b K ab = ^D a K-£ abc k bd N dc + 3K ab n b -8nGj a . (22) 

The N ab constraint equation (fT3| becomes 

2 4 
D b N ab = -D a N - e abc D b n c - -Nn a + 2N ab n b . (23) 

Finally, the spatial triad variables B k are evolved via 



D B a k = -E a b c oj„ B k - K ac B k - jB k , (24) 



and are constrained by 



2N 

e abc D b B c k = N ab B k + e abc n b B k - —B k . (25) 

IV. GAUGE CONDITIONS 

In the WEBB tetrad formalism, in addition to specifying the usual lapse function and shift vector in order to 
control the evolution of the coordinates, gauge conditions arc required to specify the evolution of the tetrad vectors. 
In our hypersurface-orthogonal implementation, the acceleration vector is given in terms of the lapse by Eq. (|8j). The 
remaining freedom is the evolution of the angular velocity (with respect to Fermi- Walker frame dragging) 3-vector 
Lo a . In this section, we derive the equation for the lapse function appropriate to our CMC hypersurface condition and 
the elliptic equations which preserve the 3D Nester gauge conditions on the spatial triad during the evolution. We 
assume for the present that the shift vector is given; discussion and derivation of the shift equations is deferred to 
Sec. EE 

A. CMC slicing 

In order to fix the lapse and the corresponding foliation of spacetime, we require that the trace of the extrinsic 
curvature be a positive constant in both space and time, 

K = const. > 0, (26) 

which with our sign convention means that the spatial hypersurfaces go to future, rather than past, null infinity. Then 
Eq. (fT8|) becomes the following linear, elliptic equation for the lapse: 



-(D a - 2n a )D a + KabKab + ^K 2 + 4ttG{ P + a aa ) 



a = 0. (27) 



We notice that the operator (D a — 2n a )D a = V a V a is equal to the covariant Laplacian on the three-geometry, where 
here and in the following, V a denotes the covariant derivative with respect to the three-geometry. It follows that the 
operator inside the square parenthesis is (formally) positive if the strong energy condition, p + <j aa > 0, holds. 

B. 3D Nester gauge 



After imposing hypersurface orthogonality, the remaining freedom in the choice of the tetrad is described by local 
rotations of the spacelike legs. Nester [25| proposed a set of natural gauge conditions in order to fix this freedom up 
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to a global rotation. These conditions are conveniently formulated using differential forms. Let 9 a , a = 1,2, 3, be the 
co-frame associated with B a ; that is, Q 1 , 9 2 , 9 s are the basis of one-forms defined by 

9 a (B b ) = S\. 

The Cartan structure equations (see for instance (ioj 1 ) imply the following relation between 9 a and the connection 
coefficients T acc i'. 

d9 a = T a cd 9 c A 9 d . 

The 3D Nester gauge condition is 2 

dE + 511 = 0, (28) 
where the one-form S and the three-form II arc defined by 

5 = i B J6 a = T\ b 9\ 

II = -6 ab 9 a A d9 b = -T bcd 9 b A 9 C A 9 d , 
2 2 

respectively. Here, 6 denotes the co-diffcrential which is defined as the formal adjoint of the exterior derivative 
operator d with respect to the inner product (v, if) = J s v A *r/ for p- forms v and rj. In terms of the Hodge dual *, it 
is explicitly given by <5 = (— l) p * d* when acting on a p-form. Provided suitable boundary conditions are specified at 
the boundary of S t , the two terms arising in Eq. (|28p are mutually orthogonal with respect to the product (•, •), and 
the 3D Nester gauge condition is equivalent to 

cG = 0, d*n = 0. (29) 

In our notation, we have S = 2n a 9 a and *n = N . Therefore, assuming that T, t is connected and that its first 
cohomology vanishes, the 3D Nester gauge implies 

n a = D Q $, N = const. (30) 

for a function $ on S t . We choose the constant value of N = 0, since this tends to minimize unnecessary twisting 
of the triad vectors from one spatial point to another, which is particularly appropriate in the context of asymptotic 
flatness. The existence and uniqueness of a set of orthonormal three-frames globally satisfying the 3D Nester gauge 
can be related to the existence and uniqueness of a nowhere vanishing spinor field satisfying the 3D Dirac equation, 
see Ref. [42j and Appendix [C] 



C. Properties of the evolution system in the CMC and 3D Nester gauge 



The gauge conditions we have adopted have the following nice properties: they fix the gauge variables K, N and 
n a (the last up to a gradient). As a consequence, the evolution system reduces to the advection equation l]24[) for the 
triad vectors B a and the Maxwell- like equations (fTol [T7)) for the symmetric, traceless parts of K ab , N ab . This system 
constitutes a quasilincar, symmetric hyperbolic evolution system, assuming that the gauge variables a b = D b (loga), 
n b = -Db<f>. and uj b arc known. The gauge variables are determined by a set of coupled elliptic equations obtained by 
requiring that the gauge conditions be preserved during the evolution. The lapse is obtained by solving the elliptic 
equation (f!?7|) . which is a consequence of setting the time derivative of K to zero, preserving CMC slicing. Next, we 
note that in the 3D Nester gauge, the Hamiltonian constraint (|2Tj) yields an elliptic equation for the scalar potential 
<& which might be used to obtain n b . Finally, the requirement that the time evolution preserve the 3D Nester gauge 
condition yields an elliptic system of equations for 2aD$<& and au> b . This system can be derived as follows. Using the 
commutation relation [<9 t — £p,D a ]& = —a(K ac + e abc uj b )D c & in Eq. (|20|) . with n a = D a &, we obtain 



D a (2aD <i>) + e abc D b {au c ) = -a 



SabcKbdNdc - K ab (3n b + a b ) + - Ka a + 8nGj a 



(31) 



2 This condition has a natural generalization to higher dimensions, see |4ll . 
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In addition, with a a = D a (logQ!) and N = 0, Eq. (fT9|) can be rewritten as 

D a (auj a ) = -aN ab K ab . (32) 

As shown in Appendix [C] Eqs. ([3T1 132|) arc equivalent to the inhomogcneous 3D Dirac equation, and they can be cast 
as a second-order elliptic system for appropriate potential fields. This elliptic formulation will be presented explicitly 
in the next section, after the conformal decomposition of the fields has been performed. 

Summarizing, in the CMC and 3D Nester gauge conditions, our system reduces naturally to a symmetric hyperbolic 
evolution problem for the fields B a , K ab , N ab which is coupled to an elliptic system which determines the gauge 
variables a, 4> and oj a . Appropriate conditions on the evolution of the spatial coordinates, and the corresponding 
equations for the shift vector, will be discussed in Sec. IV Bl The well-poscdness of the full, coupled hyperbolic-elliptic 
system depends on the structure of the elliptic equations for the gauge variables, and the choice of the shift. A well- 
posed Cauchy formulation has been obtained by Andersson and Moncrief [30j for a somewhat similar metric-based 
system. 



V. THE CONFORMALLY RESCALED EQUATIONS 

With compactified spatial coordinates and asymptotically null spatial hypersurfaces, so that .y + is at a finite 
coordinate radius the physical spacetime (coordinate) metric g^ v is singular at J? + . The spatial metric must 
generate infinite spatial distances over a finite range of coordinates, and the lapse becomes infinite as the hypersurface 
becomes null. In the conformal approach to asymptotically flat spacetimes, pioneered by Penrose [43j], these infinities 
are absorbed into a conformal factor D, such that the conformal metric g^, defined by 

is finite at with Q = at + and positive everywhere in the interior. The physical manifold M can be considered 
as part of a larger conformal manifold M containing the boundary of the physical manifold DM as a regular null 
hypersurface and in which the CMC hypersurfaces are completely spacelike. The conformal metric is taken to be 
regular at ^ + , with the precise amount of smoothness there open to some debate (see Sec. HI . 

In the 3 + 1 tetrad formalism, the above scaling of the coordinate metric under a conformal transformation corre- 
sponds to the lapse and the coordinate components of the spatial triad vectors scaling as 

a = rr 1 ^, B a l = VtB a \ 

The coordinate components of the shift vector are unchanged by the conformal transformation. We can then write 
for the directional derivatives 

D = QD , D a = ilD a . 

Note, however, that the triad components of the shift vector transform as j3 a = il(3 a . 
From the definition of the tetrad connection coefficients, Eq. (QJ, 

Tq/37 = + D a nrjj3 7 — DpQr)^. 

This translates into the following conformal transformations of the dyadic decompositions of the connection coeffi- 
cients: 

N ab = rtN ab + e abc D c n 7 K ab = ilK ab - D Q S ab , 

from which 

N ab = nN ab , N = UN, n b = nh b + D b n, (33) 
K ab = nK ab , K = QK - 3AA (34) 

and 

a b = ila b — D b il, a b = D b (log a) , uj b = £l& b - (35) 

The choice of conformal scaling for the energy-momentum tensor is somewhat arbitrary. The physical tetrad 
components of a radiation energy- momentum tensor are expected to scale as f2 4 , two powers from the inverse r 2 
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scaling of divergent radiation in a Minkowski frame, and two rcdshift-time dilation factors due to our tetrad frame 
becoming asymptotically null. Therefore, we will define the conformal energy-momentum tensor by 



T a /3 



P = T"00, 



Jb 



Cab = Tab- 



Ill principle, there is quite a bit of freedom in the choice of the conformal factor O. Any function which is positive 
over the physical manifold and vanishes in a suitable way at would do, since a change in Q can be compensated 
for by a change in the determinant of the conformal metric. One possible choice would be to demand that f2~ 6 equal 
the determinant of the physical spatial metric, so the determinant of the conformal spatial metric is 1. However, 
such a choice would be highly coordinate dependent. What seems to us to be a much better choice, given our use of 
the 3D Nester gauge, is il = e*, where $ is the Nester potential for the physical space, see Eq. (|30|) . This has the 
great simplifying consequence that n a = [26j . and together with N = this means that the conformal connection 

coefficients N ab = N ab are symmetric and traceless. Another consequence is that the triad vector fields, as vector 
fields in the conformal 3-space, are divergenceless: 



(36) 



V k B a k = -2h a = 0. 
The covariant conformal derivative of an arbitrary vector field v = v a Ti a is 

X? a Vb = D a v b - N ac v d e cdb , 

and V a t> a = D a v a - Therefore, in the conformal Laplacian of a scalar or the conformal divergence of a vector field, the 
covariant derivatives can be replaced by directional derivatives. Furthermore, the conformal Ricci tensor is 



R ab = -D c N d(a e b)cd + 2N ac N bc - S ab N cd N cd . 



(37) 



The Ricci scalar has the simple form R = —N ab N ab and is negative definite. 

The Hamiltonian constraint equation (|2f [) then becomes a relatively simple elliptic equation for fi: 



n D c D c n 



D c n) (D c n 



<r 2 



K cd K cd + N cd N cd ) + 4irGn 4 p 



(38) 



With the Nester-gauge inspired choice of 0, and with the time derivatives of f2 expressed in terms of K, the 
evolution equations for the conformal quantities are 



D n 

D B a k 

D N ab + D c K d ( a e b ) cd 
D K ab — D c N d ( a e b y d 



K - QK 



-K ab B b k 



£acd^c.B d k —B a k , 



K A 



K - 



2K c[a N b)c - Y Nab + 2e cd(aN b ) c U>d + -r 

-2N ac N bc — —R'ab + 2s cd ( a K b ) c cb d 



£cd( a K b ) c D d a 



D 



(a (a W 6 ) J 



TF 



(39) 
(40) 

(41) 



!^ « - 2 
-V a V b a- - 
a 11 



v a v b n + —k ab 



87rGf2 2 <7 a (, 



TF 



This system is subject to the constraint equations 



~-abcD b B c k 



N ab B b k = 0, 
D b N ab = 0, 



D b K ab + e abc K bd N dc 



^(D b fl)K ab 



-8irGn 2 j a . 



(42) 

(43) 
(44) 
(45) 



The conformal lapse a is determined by transforming the CMC hypersurface condition, Eq. (|27p . to conformal 
variables and then using the Hamiltonian constraint, Eq. (|38[) . to eliminate the most singular terms at , giving 



flD a D a & - 3(D a Q)D a & + (D a D a Q)& - - (^N ab N ab + 3K ab K ab ^ a = 4^Gfi 3 (3p + a cc ) a. 



(46) 
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We see that in principle, ft can be found in two ways: by solving the elliptic equation (|38|) of the Hamiltonian 
constraint or by integrating the evolution equation ([59]) . The elliptic equation is degenerate at ^ + due to the 
hypcrsurface becoming null. The apparently singular terms constrain the asymptotic behavior of the solution, but in 
practice cause no numerical problems fLU [TtJ ■ Likewise, we do not expect the apparently singular terms in the lapse 
equation (|46[) to impede the use of standard elliptic solvers. The time evolution equation for ft could prove useful in 
reducing the frequency of solving its elliptic equation. However, as yet we have no equations for K and w a . Elliptic 
equations for these quantities come from requiring preservation of the Nester gauge during the evolution, as shown in 
the following subsection. 



A. Determination of K and ui a 

Now we derive the elliptic equations for the gauge fields K and ui a . These equations arc based on Eqs. (f3"2")l and ([5T|) 
from Sec. II V Cl which yield, in terms of the rescaled quantities, 



D a I \aK \ + e a bcD b {aCj c ) = afl a , (48) 



D a (aQ a ) = aN ab K ab: (47) 

-I- F„u. 

where 

fra = K ab D b (log a + 2 log ft) - e abc K bd N dc - 87rGft 2 j a . (49) 

While Eq. (|47l) depends only on the rescaled quantities, jl a in Eq. (|48[) contains the apparently singular term involving 
the gradient of log ft. However, jl a is, in fact, finite at ,y + . Indeed, using the rescaled momentum constraint, Eq. (|45l) . 
we may rewrite 

ajjb a = b b (ak ab ) (50) 

which does not involve the conformal factor. However, we prefer Eq. (|4"§)). since it does not contain derivatives of 

dynamical fields. In Sec. IVI1 we derive an explicit regular expression for the term K ab D b log ft at . In Appendix [Al 
we show that second radial derivatives of K ab are logarithmically singular at J? + whenever outgoing radiation is present 
there. 

In order to cast Eqs. (|471 I48p into a second-order elliptic system, we make the following ansatz which is motivated 
by the considerations at the end of Appendix [Cl 

2 - 

-aK = -D a ( a , auj a = D a ip + e abc D b C, c , 

where ip and ( a are a scalar function and a vector field, respectively. Plugging this into Eqs. (|T71 and using the 
commutation relation [D a ,D b ] = s ab dN c dD c , we obtain the elliptic system 

-b a b a i> - N ab b a c b = &N ab K ab , (51) 

-D b D b ( a + N a bb b ip - e abc N cd b d ( b = aft a . (52) 

When supplemented with appropriate boundary conditions at J^ + , this system yields a unique solution for {ip,Ca)- 
We choose homogeneous Dirichlet boundary conditions for ip and the tangential components of Q a . The normal 
component of C, a , however, is determined by the boundary condition from Eq. (|A31[) which guarantees that the cross 
sections ^ + remain metric two-spheres. This translates into a mixed boundary condition for the normal component 
of Ca (sec Sec. |VH] below). 



B. Shift Condition 



The shift governs the evolution of the spatial coordinates. Therefore, keeping the coordinate location of , where 
ft = 0, at some fixed coordinate radius R = R+ imposes a boundary condition on the shift. The equation for the 
evolution of the conformal factor is 

ab Q n = (d t - p k d k ) n = -a(j- ^k) , (53) 
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with K finite at J? + . The requirement that the conformal factor remain zero at J? + , <9tf2 = 0, implies that /3 fc dfcf2 = 
(3 a D a £l = Ka/3, which gives a constraint on the normal component of the shift. The outward normal vector field s a 
with respect to the conformal geometry is proportional 3 to — D a Q,, and the normalization factor can be determined 
by evaluating the Hamiltonian constraint, Eq. ([55)1 . at J^ + . This gives 



K ' 



(54) 



so that the normal component of the shift at J" + is s a j3 a = —a. It also makes sense to require that the components 
of the shift tangent to J? + vanish. This means the coordinates on J? + are propagated along the null generators of 
J^ + . The complete shift boundary condition is then 



A, 



K K 



(55) 



Satisfying this boundary condition is most straightforward if it is imposed as a boundary condition on a vector elliptic 
equation for the shift. This spatial coordinate gauge fixing may need to be done carefully, and here we mention a few 
possibilities and their main features. 

The first possibility, the spatial harmonic gauge, is motivated by the metric-based formulation proved well posed by 
Andersson and Moncrief [30(. With CMC foliations, this gauge casts the York version 44| of the standard Arnowitt- 
Deser-Misner (ADM) evolution equations [45| into a nonlinear system of wave equations which is coupled to elliptic 
equations for lapse and shift. Besides bringing the evolution equations into hyperbolic form, the spatial harmonic 
gauge has additional nice properties which were noted in [30j and exploited in their proof. Namely, the terms involving 
derivatives of the three-metric in the elliptic equation for the lapse are eliminated. Similarly, the terms involving second 
derivatives of the three-metric in the vector elliptic equation for the shift cancel out. These properties were important 
for consistently solving the elliptic equations to agiven degree of smoothness in the rigorous proof of well-posedness 
for the mixed hyperbolic-elliptic system in Ref. |3(| ■ While their analysis was carried out on compact CMC slices in 
the absence of conformal rescaling, it motivated the choice of the spatial harmonic gauge for the conformal spatial 
geometry in the conformally compactified scheme of Moncrief and Rinne [l3| . 

In terms of the conformal metric h^ the spatial harmonic gauge sets 



V 



k — 



h' J 



1 ij J- i 



(56) 



where f fe y- and T k ij are, respectively, the Christoffel symbols of the conformal three metric hij and of a given, possibly 
time-dependent, reference metric hij on S t . Requiring this gauge condition to remain satisfied during the evolution 
implies an elliptic equation for the shift, derived from dtV h = 0. 

In general, there is no simple way to express the Christoffel symbols in terms of the triad connection coefficients. 
However, in the conformal Nester gauge, there is a very simple expression for the trace T fc = h %3 T ij in terms of the 
conformal triad vectors. A straightforward calculation gives 



-dih 



D a B a k - (ViBj)B a 



In the conformal Nester gauge, V,i? a 4 = (see Eq. (|56")l ). so 



f k = -D a B a k . 



(57) 



The elliptic equation for the shift comes from substituting Eq. (|40|) into 



-D„ 



(d t ~ ^d t )B a k - {DaP^B^ 



(d t - P%)Bj] (djB a k + 2/>'„'l ' ,.) - B a l B a i(d t - fi l d t )T k 



0. 



3 The minus sign in the expression —D a Q comes from the fact that Q is strictly positive in the interior of St, hence, its gradient points 
inward from J^+. 
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After using Eqs. P51 - |5U| to eliminate derivatives of dw a and aK ab , the result is 



D a D a /3 k + 2B^T k tJ D a f3 l 



tafia - ^D a (aK] 



B„ 



+ 2aK ab (b a B b k + ll/lh'Y'-,) - B a l B a i(d t - pd t )T l 



0. 



(58) 



It may be significant that the spatial harmonic gauge eliminates derivatives of the triad vectors from the triad 
"Laplacian" operator D a D ai which is the true scalar Laplacian operator in the conformal Nester gauge. In the spatial 
harmonic gauge, 



D a D a = B a *BJ 



2 



D a B a k ) = BaBj 



dx l dx : > 3 dx k 



In all of our elliptic equations except the shift equation, this completely eliminates derivatives of the conformal triad 
vectors. Notice also that the elliptic equation for the shift required for preserving the spatial harmonic gauge, Eq. (|58[). 
does not contain second or higher derivatives of B a k . An appropriate choice for the reference metric of the spatial 
harmonic gauge may well be the trivial one h%j = Sij, for which the Christoffel symbols vanish. This can accommodate 
multiple black holes, as shown in the conformally flat initial data calculations of Ref. [l7| . 

We now consider alternative coordinate gauges, based on the motivation of choosing the shift such that the spatial 
metric components become time independent as the physical and presumably the conformal spacctimc become sta- 
tionary. The spatial harmonic gauge may well have this property, but a more direct way of achieving the goal is to 
choose the shift to minimize a functional quadratic in time derivatives of the conformal metric. The two options of 
this type presented here are contingent on establishing analytically, or by successful numerical trials, that they do not 
lead to instabilities. We focus on minimizing the time dependence of the conformal metric, rather than the physical 
metric, because our conformal factor is a coordinate scalar determined by a coordinate scalar equation. 

The time dependence of our conformal metric is determined by, though not equivalent to, the evolution equation (|40)) 
for the conformal triad fields, which is regular at J? + . Using the constraint equation we can rewrite Eq. (^0)) as 



d t B a l = - 



DaPb + £acdPcN db + a (R 



B h 



(59) 



The triad components of the shift depend on the conformal rescaling, hence the j3 a . Eq. (|59"|) should supersede Eq. ([^0]) 
as the evolution equation for the triad vectors when the shift is calculated as the f3 a rather than the /3\ Contracting 
with the dual one-forms 6 b i gives 



-D a l3 b - e ac d PcNdb ~ Ot ( K ab + E acb U), 



F ba — i d t Ba 

One should bear in mind that the determinant of the conformal metric is not time-independent in general for the 
Nestcr-based choice of conformal factor. Therefore, the straightforward way to minimize the time dependence of the 

conformal metric is the "minimal conformal strain" condition minimizing J[(dth^)(dth ki ) hikhji/w \fh] d?x, where w 
is a weight function. Since h %3 = B a l B a J , this is equivalent to minimizing the functional 

(60) 



h\Pa] = J F(ab)F(ab)WTj, 

where i) = \/h,d 3 x is the natural volume element of the conformal geometry. The trivial choice of weight function, 
w = 1, seems as good as any. The symmetrization of F ab in Eq. (|60p is what makes it essentially equivalent to the 
conventional minimal strain condition in terms of h lJ . 
With 

F (ab) = -V( ^b) - aK ab , 
where V a is the conformal covariant derivative operator, minimizing Ii gives the elliptic equation 

- V b F (ab) = - (V 6 V 6 /3 a + V 6 V a ^ b ) + V b {aK ab ) = 0, (61) 
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and the boundary variation 5/3 a F^ a ^Sb vanishing is consistent with the Dirichlet boundary condition (|55p at 
discussed above. As in the previous subsection, we prefer avoiding radial derivatives of the fields in the source terms 
of the elliptic equations. For this reason, we use the momentum constraint equation (j45|) . and substitute 

\7b(aK a b) = & 

As mentioned before, a regular expression for the term K a bF>b log ft at will be derived in Sec. I VII We note that 
expanding the covariant derivatives gives 



V fc V 6 /3 Q = D b D b Pa + 2e acd N b dDbl3 c + N ab N bc /3 c - N bc NbcPa 

and 

V h V a /3 h = V a (V 6 &) + Rabh = Da (pbPb) + RabPb, 

where the constraint Eq. (|44j) has been used and R ab is given by Eq. (|37|) . 

Smarr and York [46j argue that a different shift condition, the "minimal distortion" condition, most cleanly sup- 
presses that part of the time dependence of the spatial metric associated with diffeomorphisms of the spatial coor- 
dinates. This minimizes the time dependence of what they call the conformal spatial metric, = (/1 1 / 3 ) /i 4 -? , by 
minimizing the action J[(dth l ^)(dth ki )hikhjg\^h] d 3 x. Their rationale for minimal distortion rather than minimal 
strain is based on using the Hamiltonian constraint to determine the determinant of the spatial metric, but we make 
a quite different kind of conformal decomposition. For completeness, we note that the minimal distortion shift can be 

implemented in our formalism by defining B a l = B a t /(detB) = B a l /(det B) and F a b = F a b — \& a bF cc . The minimal 
distortion version of the action is 

Wa] = J F {ab) F {ab) fi. (62) 

The differential equation for the shift is then 

= -V b F {ab) = i (VbVbpa + V 6 V a /3 6 ) - iv a V 6 /3 6 + Vfc (&K ab ) 

= \ {VbVbPa + iv a V b /3 6 + RabPb^J + V b (aKab) , (63) 

and the boundary variation is 5(3 a F^ a b)Sb, which vanishes for the Dirichlet boundary condition ()55[) . 

All of the above forms of the shift equation are regular at ,f + and should give rise to a well-posed elliptic boundary 
value problem on a CMC hypcrsurfacc (with other variables fixed) as long as the only boundary is at . If black 
holes are present, there should also be inner boundaries at excision surfaces inside or at the apparent horizons, but 
we defer discussion of inner boundary conditions to a future paper. 



K ab Db(loga + 2 log ft) - 8nGn 2 j a 



-D a ( aK 



VI. REGULARITY CONDITIONS AT ,f + 

The evolution equations (|39l HOl BT]) for ft, B a k and N a b are completely regular at J^ + . However, Eq. (|4"2"T) for the 
evolution of K a b has the apparently singular term 

1 /- ~ K ~ 

s ab = ^ lv a v b n + —K ab 

at J? + . In this section, we derive the appropriate regularity conditions which ensure that S a b is finite at To 
analyze the limit of the covariant Hessian S a b = V a Vbft = V( a Vb)ft at J^ + , we introduce the tensor ^f a b — S a b — s a Sb 
projecting tangent to ^0+ and the 2D extrinsic curvature of 

Kab = laclbd^cSd = Jac^cSb, (64) 



± 1 
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with a similar sign convention to that of the 3D extrinsic curvature, such that positive for expansion of the 

outward normals. With s& = — (3/if)V&f2 and smoothly extending as the tangent vector to the normal spacelike 
geodesic in a neighborhood of J r +. decomposing E a b into components along and perpendicular to s a gives 

S q6 = SaSbScSd^c^d^ + S a S c 76dV d V c Q + S b Sdlac V c V d Sl + 7ac7bdV c V d f2 

= s a s 6 s c V c (^dVd&j - yK a 6, (65) 
since s c V c Sd = 0, k oc s c = and £1 = 0. Take the trace of Eq. (p5|) to get 

V c V c f2 = ,S C V C (sdVd^J - y/C, 

and the gradient of Eq. (|3"5j) gives 

§b^ab = ^S a DbDb^ = gS a S c V c (s d V d^j - yKSa, 

from which 

S c V c (s d V d fi) = S c Sd^ab = -UK/6 

and 

Sab = -y {kab + ^KS a S b ^ . (66) 

Therefore, we find for the singular term in Eq. (|42[) 

QS a b = (S ao ) TF + y^ ab = — y — ~jlabk — ^ ab 
To ensure that S a b is finite, we see that the condition 

kab = k ab - - J a bk = K a b (67) 

must hold at J" + . and both 

K ab h = 0(Q) (68) 

and 

kab - laclbdKcd = 0(ft) (69) 

at J+. 

The first of these conditions, that only the transverse-traceless part of K a b can be non-zero at ^ + , is what is 
required to resolve the apparently singular terms in the sources for the elliptic equations discussed in Sec. [Vj As 

we show at the end of this section, the limit of Q~ 1 K a b§b at J'^ can be derived from the momentum constraint 
equation j27| . The second condition, Eq. (|69p , is known as the zero-shear condition, and physically is a constraint on 
the initial data that limits the amplitude of incoming waves in the neighborhood of ^ + and ultimately at past null 
infinity. There is an extensive literature (see [27]]) showing how these conditions relate to the existence of ^ + as a 
regular null hypcrsurface in the conformal spacctime. 

While Eq. (|66p can be inverted to get an expression for k a b at in terms of directional derivatives of ft, it is 
better to evaluate k a b from its defining expression 

kab = 7ocV c S 6 = J ac (^D c S b + £bde.SdN c ^j , (70) 

which does not contain any second normal derivatives of Q. Since the directional derivative of s b is tangent to J^ + , 
k a b at is obtained by substituting §b = —(3/K)DbCl into Eq. ((70)) . The trace of the 2D extrinsic curvature is then 

3 - ~ 

k = JabDbSa = -—JabD b D a n. (71) 
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Away from J? + , s b is found by integrating the geodesic equation. 

In the following, we compute an explicitly finite expression at J'^ for the apparently singular term S a b when the 
regularity conditions (|68f and (|69|) are satisfied. The direct approach using l'Hopital's rule would require going to 
higher order in the asymptotic expansions, which is awkward in the context of the tetrad equations. The Nester 
gauge conditions, as elliptic equations which are regular at J? + , impose no local constraints on the behavior of the 

N a b near or at . On the other hand, the asymptotic expansions are rather straightforward when dealing with 
coordinate components in certain special coordinate systems. We derive in Appendix [A] the asymptotic expansion 
of the conformal factor in Gaussian normal coordinates, and from this, expressions for the coordinate components 
Sij. Since S is a 3-tcnsor, it is then trivial to obtain the corresponding expression for the triad components, though 
evaluating this expression is not so easy, since it requires knowing s a away from 

An alternative approach is to impose the Penrose regularity condition that the Weyl tensor vanish at J? + . In 
Appendix [A] we use the asymptotic expansions in Gaussian normal coordinates to evaluate the Weyl tensor and 
determine the additional constraint on the asymptotic behavior of the extrinsic curvature beyond the zero-shear 
condition required to satisfy the Penrose condition. While it is still somewhat controversial, we believe that the 
Penrose condition is reasonable, in that it is only a slightly stronger restriction on the amplitude of incoming waves 
near J? + in the initial data than the zero-shear condition, and is preserved in the subsequent evolution. 

A standard expression for the electric part of the physical Weyl tensor, in which time derivatives have been eliminated 
using the evolution equations (as in Ref. [13j . Eq. (5.1)), is 



E a b = CoaOb = [Rab - K ac K cb + KK ab — \ixGa ai 



; s K ~ 

Rab — KacK c b + ~^K ab 



4irGcrai 



'TF 



(72) 



The conformal invariance of the Weyl tensor means that the conformal Weyl tensor, i.e., the Weyl tensor calculated 
from the full space-time conformal metric, is E ab = Q~ 2 E ab . It is the coordinate components E^ which are truly 
invariant, and the scaling relation between physical and conformal triad vectors introduces the ft -2 factor. In our 
conformal variables, and taking advantage of the simplifications of the Nester triad gauge, we have 

-2] 



E a b — ft 
1 

n 



'E a b 

K ~ 

V a V 6 ft + — K ab 



TF 



Rab - K ac K bc - 4irGfl 2 & at 



TF 



S ab - D c Nd{ a £b)cd + 



2N ac N bc - K ac K bc - 4^Gft 2 a at 



TF 



Note that the expression in Eq. (|72[) is not conformally invariant, since the Einstein equations used to transform the 
electric tensor to that form are not conformally invariant. The Penrose regularity condition implies that E ab = 0, 
such that 

- ~ i ; iTF 



S a b 



D r N, 



c^d{a £ b)cd 



2N ac N bc 



' 6c 



Introduced into the evolution equation (|42p for K ab , the result is 



E>oK ab = -D c N d ( a e b ) cd 



2 N ac N bc - K ac K t 



be 



K ~ ~ 1 - - 

-K ab + 2e cd ^ a K b - )c uj d + — V Q Vf,a 



-i TF 



(73) 



(74) 



3 a 

The conformal momentum constraint, Eq. (|45p . can be used to obtain more explicit information on the asymptotic 
behavior of the extrinsic curvature. Start from the equation in the form ^ -1 if a b^ — Q~ 1 (VbtyK ab = —8irGfl 2 ja- 

Decompose K ab into longitudinal and transverse parts relative to the normal vector s a at , 

Kab = S a S b (scSdKcd) + 2S( a {^b)dS c K cd j + (jaclbdK c dj ■ 

Substituting into the momentum constraint and using the properties of s a gives 

s c s d K cd \ a ^\ / s c ^bdK cd 



^SfcVb + Clk — s b D b Q 

^fisbVb + fi«; — s b D b fl 
-8nGn 2 j a . 



ScladK cd 
ft 



(ftV b - L> b ft) 



ft 



Kab S c % d K, 



cd 



> £ b n 
' ft 



I laclbdK cd 
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The purely or partially longitudinal parts of K a b/Q are finite at J^+, though we shall sec in Appendix [Al that their 
radial derivatives are logarithmically divergent. Despite this, those derivative terms vanish at since they contain 
an fi factor in front. Projecting along s a and evaluating at J? + yields 

Cl~ l S c SdK c d = Jf^ab (laclbdKcd^J = J^K cd K cd . (75) 

Since ^ c dK c d = —s c §dK c d, which is zero at J? + , K a b can be replaced by its traceless part k a b = H a b — \lab^ m Eq. (|75l) . 
The projection perpendicular to s a gives at / + 

Q~ 1 S c JadK c d = --^7ae7&dVb (jceKcd^j ■ 

Combining these results, the Q~ 1 st,Kt, a term that appears in various equations can be explicitly evaluated as 

£l~ 1 s b K ba = -— 7 ae 7hrfVf, (jecKcd^j + j^s a k bc Kbc- 

Based on the asymptotic expansions in Gaussian normal coordinates, we show in Appendix [A] that the regularity 
conditions K ab Sb = and "f a c"fbdK cd = h a b are preserved by the time evolution. 

VII. SUMMARY 

The mixed hyperbolic-elliptic system proposed in this paper for solving the Einstein equations on CMC hypersur- 
faces in a conformally compactified asymptotically flat spacetime is, we believe, worth testing as a way of improving 
the accuracy of numerical calculations of gravitational radiation emitted in mergers of black holes. There is no need to 
use artificial boundary conditions at a finite radius when the computational domain extends all the way to future null 
infinity, and asymptotic gravitational wave amplitudes are calculated directly, without the need for any extrapolation 
procedure. A purely hyperbolic system would possibly be more computationally efficient, since there are no elliptic 
equations to solve, but requires extra care in constructing a numerical scheme to suppress errors in the additional 
constraints among the additional variables such a system entails. 

Our representation of tensors by tetrad, rather than coordinate, components is perhaps more controversial. We fix 
the timelike vector of the tetrad by requiring it to be orthogonal to the CMC hypersurfaces of constant coordinate 
time, but there is gauge freedom in the orientation of the triad of spatial vectors within the hypersurface, in addition to 
the gauge freedom in the evolution of the spatial coordinates. A key part of our scheme is combining what wc consider 
to be the "natural" choice of triad gauge, the 3D Nester gauge, with the choice of conformal gauge (definition of the 
conformal factor), such that the non-trivial spatial triad connection coefficients for the conformal geometry of the CMC 

hypersurfaces are equivalent to the five independent components of the traceless symmetric N ab - The hyperbolic part 
of our evolution system for the conformally rescaled variables consists of the coupled evolution Eqs. (|41[ for the 

N a b and the traceless part of the conformal extrinsic curvature K a b, respectively, plus the simple advection Eq. P0|) 
(or ([59]) when f3 a is used) for the nine coordinate components of the conformal triad vectors, B a l . This is a total of 
nineteen evolution equations. 

The remaining six tetrad connection coefficients - the triad components of the acceleration and angular velocity 
(relative to Fermi- Walker transport) of the tetrad frames - are determined by elliptic equations: Eq. (fH)]) for the 
conformal lapse, whose logarithmic gradient is the conformal acceleration, and Eqs. (|5 1 [) and ([5^)1 for a scalar potential 
ip and a vector potential C Q , which together determine the conformal angular velocity uj a and the trace of the conformal 
extrinsic curvature K. The conformal gauge is enforced by the Hamiltonian constraint Eq. (|38[) for the conformal 
factor fl, but some interpolation may be possible using the evolution Eq. (j3"5)) . The coordinate gauge conditions are 
fixed by elliptic equations for the triad components of the conformal shift. We suggest three options, Eqs. (J5HJ), (foTj) . 
and (|63|) . but some experimentation will be necessary to determine which of these, or perhaps some quite different 
scheme, works best in practice. We do not propose explicitly imposing the momentum constraints during the evolution, 



since this does not seem to be necessary for stability [14j, and we similarly allow the DbN a b = constraint to evolve 
freely. We analyze the propagation of errors in these and other constraints in Appendix [B] 

No boundary conditions are necessary at J'^ for the evolution equations, which is a major advantage of the CMC 
hypersurface condition. Boundary conditions are required for the elliptic equations. Fundamental to the conformal 
decomposition is the boundary condition == on the conformal factor. As discussed in Sec. IV Bl the boundary 
condition at ^+ on the shift vector is (3 a = —as a , where s a = —(3/K)D a fl is the unit normal to in the 
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conformal 3-space. This ensures that J> + is at a fixed coordinate radius R + and that the angular coordinates in ,f + 
are propagated along the null generators of These angular coordinates are related to the nominally Cartesian 
computational coordinates in the same way as in flat space. 

The boundary conditions for the Nester gauge potentials ijj and £ a introduced in Sec. IV Al are naturally chosen to 
be homogeneous Dirichlet conditions for ip and the part of ( a tangent to : 

^ = 0, % b (b = 0, (76) 

since we do not want the conformal triad angular velocity Q a to be unnecessarily non-zero. For a flat conformal spatial 

geometry, N a b and therefore the source terms for ip are zero, which should make "0 — globally. Furthermore, in a 
spherically symmetric spacetime, the boundary conditions should allow £ a to be spherically symmetric, which means 
it does not contribute to U3 a . 

In order to keep the intrinsic geometry of the ^+ 2-surface that of a true sphere, with an inverse circumferential 
radius £ constant in time, as discussed in Appendix [5] the boundary condition (|A31|) on K must be satisfied. In 
terms of the Nester vector potential, this becomes Z? a Ca = —ak, or 

SbD b (5 a ( a ) - ks a C a = -OiK. (77) 

While we use Gaussian normal coordinates to derive this condition, it is valid regardless of the computational coordi- 
nate system since it is expressed in terms of spatial coordinate scalars. The area of the J* + 2-sphcre is then constant 
in time, which means that the expansion of vanishes. Note that £o is not necessarily equal to the inverse of the 
coordinate radius of R+. The geometrical definition of £o is that the 2D scalar curvature of ,f + is 2£ 2 . In 
Sec. I Villi we use the zero expansion result, plus the fact that our boundary conditions on the shift maintain a simple 
relation between our computational coordinates and standard polar angle coordinates on to obtain a simple 
expression for the Bondi news function in terms of our variables. 

Finally, we consider the boundary condition on the conformal lapse. We can make our time coordinate correspond 
to retarded Minkowski time at ^ + by imposing the boundary condition 

a = — = a (78) 

on Eq. (pQj)) for the conformal lapse, since CMC slicing of Minkowski space has a ~ Kr/3 as r — > oo. Since £ is 
uniform on all of the hypersurface to the future of the initial hypcrsurface, so is do- 

The question of what boundary conditions to impose on an excision surface on or inside the apparent horizon of a 
black hole is left open for the present. 



VIII. DISCUSSION 



An important issue when formulating the equations on CMC hypersurfaces is dealing with singularities, real or 
apparent, which are associated with the hypersurfaces becoming asymptotically null as they approach The 
elliptic equations for the conformal factor and the conformal lapse a are genuinely singular there. However, the 
singular terms in these equations do not cause any singularity in the solutions; instead they force any solutions which 
satisfy the boundary conditions at to have very particular asymptotic behaviors. These asymptotic behaviors, 
additionally constrained by the usual "zero-shear" and Penrose regularity conditions, were discussed in Sec. I VII We 
also found it useful to explore the asymptotic behavior at greater depth in a special coordinate system, Gaussian 
normal coordinates (Appendix |A"]) . These coordinates are not suitable as global coordinates in an actual numerical 
calculation, but their simplicity in a neighborhood of facilitates drawing general conclusions about the asymptotic 
behavior of scalars and complete tensors. In this way, we are able to analyze the apparently singular terms that appear 
in some of the evolution equations, show that they arc in fact finite at J? + , and evaluate their limits at J> + . However, 
the conformal factor, the conformal lapse, and the conformal extrinsic curvature are not completely smooth at 
on a CMC hypersurface when any outgoing radiation is present. The fourth normal derivative of f2 and the second 
normal derivative of the conformal extrinsic curvature are logarithmically infinite there. This is a gauge artifact of the 
CMC slicing condition, has nothing to do with our spatial gauge conditions, and is not prevented by any physically 
allowable regularity conditions. In principle, it should slow the convergence of a pseudo-spectral numerical scheme 
or a high-order finite difference scheme, but the impact may not be severe. The only apparent way of avoiding this 
potential problem is to make the lapse evolve dynamically as part of the hyperbolic system. 

It is instructive to give a more detailed comparison of our system with the coordinate-based system of Moncrief and 
Rinnc il3j. They have a set of eleven evolution equations for the six components of the conformal spatial metric and 
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the five components of the traceless part of the canonical momentum ir trl K The traceless canonical momentum is, 

within a factor of the square root of the determinant of the conformal metric, equivalent to our K a b- The Moncrief- 
Rinnc evolution equations are second order in spatial derivatives of the metric and much more complicated than 
our evolution equations, but there are fewer of them, both because our system is first-order and because the spatial 
triad vectors have three more degrees of freedom than the conformal metric. Since Moncrief and Rinne also impose 
the CMC hypersurface condition, their equation for the lapse is essentially the same as ours, and there is the same 
breakdown in smoothness due to polyhomogencous terms. Their conformal factor is defined differently by demanding 
that the conformal scalar curvature be uniform in both space and time, but it satisfies a similar elliptic equation. 
The difference in conformal gauge does not affect the leading behavior of the conformal factor at and the basic 
structure of the apparently singular terms in the equations is the same. Since the triad angular velocity plays no role 
in a coordinate-based system, Moncrief and Rinne have three fewer elliptic equations than we do. Their variable Y 
defined in their Eq. (2.11) is equivalent to our aK, but is determined by a quite different elliptic equation. They note 
the advantages of imposing a boundary condition on Y equivalent to our Eq. (|A31[) . Their elliptic equation for the 
shift preserves the spatial harmonic condition on the spatial coordinates. 

The tradeoff between complicated spatially second-order evolution equations and more, but simpler, first-order 
evolution equations will probably make little difference in computational efficiency; however, high accuracy may be 
easier to achieve with our first-order equations. The more serious numerical efficiency issue in our formalism is the 
three extra elliptic equations. All of our variables except for the triad vectors and the shift vector are spatial coordinate 
scalars, which with the globally optimized triad orientations provided by the Nester gauge, may be better behaved 
in complicated dynamic scenarios than the coordinate components of the metric tensor and the conformal extrinsic 
curvature. In both approaches, it is critical to use highly optimized elliptic solvers. 

Our asymptotic wave amplitude is h a b = K a b 7 whose coordinate components are Xab in the asymptotic expansion 
of the spatial metric in Gaussian normal coordinates developed in Appendix [XJ To see how this relates to the Bondi 
news, we solve asymptotically for the coordinate transformation between our Gaussian normal coordinates defined on 
CMC hypersurfaces and Bondi coordinates. The inverse metric in Bondi coordinates u, x, x A , with u the retarded 
time, x the inverse of the Bondi radius r, and x A ' the angular coordinates, is 

g uu =g UA '=0, g ux =x 2 e~ 2 P, g xx = x A We~ 2p , g XA ' = x 2 U A ' e~ 2 ^ g A ' B ' = x 2 h A ' B ' , 

with the determinant of h A B equal to (sin6>)~ 2 . For a certain class of solutions of the Einstein equations, the Bondi- 
Sachs metric functions j3, W, U A ' , and h A ' B ' have expansions in powers of x away from future null infinity at x = 0, 
provided that the Penrose regularity condition is satisfied (see |29j). 

For the special class of inertial Bondi coordinates (see Sec. 6 of Ref. 0] and Ref. @), the angular metric h A > n i has 
the expansion 

h A > B > = h A 'n' + Xa'-b'X + 0(.t 3 ), 

where as in Appendix lAl h A > B i is the unit sphere metric, Xa'b' is a traceless, symmetric tensor on the unit sphere, 
and the absence of an 0(x 2 ) term is equivalent to the Penrose regularity condition in this context. Also, j3 = 
and U A ' = 0. The Bondi mass aspect M, whose average over solid angle is the Bondi energy, is defined by W = 
1 — 2M(u,x A )x + 0(x 2 ). The Bondi news is basically the time derivative of Xa'b', often represented as a complex 
scalar c whose real part is yfe 1 = ~X P ip > an d imaginary part is x^g' sm #' = X° v ' I sin.0'. The square of the Bondi news 
at a given physicalpoint on is invariant under transformations of the Bondi coordinates from one inertial Bondi 
frame to another j4[. 

In Appendix [XJ we discuss the asymptotic conformal spatial metric in Gaussian normal coordinates z, x A and 
obtain asymptotic expressions for the lapse and shift which preserve the CMC hypersurface and Gaussian normal 
coordinate conditions. While these spatial coordinates are not the computational coordinates, the shift satisfies the 
same boundary condition at . The boundary condition on the conformal lapse, as discussed in Sec. IVIlT ensures 
that the coordinate time t on the CMC hypersurfaces asymptotically corresponds to Minkowski retarded time. The 
coordinate transformation from these coordinates to the Bondi-Sachs coordinates can be constructed as a power series 
in the distance z from starting from u = t and x A = x A , by requiring that the transformed metric have the 
Bondi-Sachs form, with g uu = g UA ' = 0. The result is 

u = t+^ + ^z 2 + 0(z 3 ), x A ' = x A + 0(z 3 ). 
2ao oao 



Then 
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and equating the determinants of the two sides gives 

x = W(l + 0(z 3 )) = ^z + 0(z 2 ). 



Note that dx/dt = 0(z 2 ), since our boundary conditions keep d£/dt = 0. From the expansions of h A ' B ' and h c ° 
(Eq. IA5[) , we see that 



6 . . 

Xa'b' = -t^Xab- (79) 



The expansion of .y + clearly vanishes. From 



or oz at 



= O(z), and from 

U A ' = 0(z 2 ). Therefore, all of the conditions for an inertial Bondi frame are satisfied, and the Bondi-Sachs "news" 
is Xa'b' = 9 U Xa'b' ■ The vacuum Einstein equations in Bondi-Sachs coordinates give, since d u = dt, 

- _I/,A'C',>0- , ± , I^A'B' 
Qy, ~ ^ n XA'B'XC'D' + 4 X |A'B' 



(80) 



with T denoting covariant derivatives on the unit sphere |47| . The monopole and dipole parts of the term linear in x AB 
as integrated over the unit sphere vanish identically, so this term contributes to neither the Bondi energy nor the Bondi 
momentum. Eq. (|80[) defines the precise physical meaning of Xab as the 2-tensor describing the asymptotic amplitude 

of the gravitational waves. This 2-tensor is equivalent to the 3-tensor whose triad components are h a b = K a b, a s 
projected into =/'+. 

Our extraction of the Bondi news is actually quite a bit simpler than it is in Cauchy-characteristic matching or 
Cauchy-characteristic extraction formalisms Q. In the characteristic methods, the matching or extraction is done on 
a world tube in the interior, and the constraints controlling the conformal Bondi frame are integrated outward on the 
null hypersurface, generically resulting in a non- inertial Bondi frame at J^ + . The subsequent transformation from 
the computational frame to the inertial Bondi frame requires considerable numerical effort @ , but a more direct way 
of calculating the Bondi news from an arbitrary characteristic frame is presented and tested in Ref . Q . 

Although we contend that the approach we have presented holds much promise for high accuracy calculations 
of gravitational waveforms from 3D numerical simulations of the Einstein equations, it remains to be seen whether 
or not a stable and efficient numerical implementation can be achieved. A good place to begin is to fine tune the 
numerical solution of the elliptic equations of our method in the context of the initial value problem and then evaluate 
whether or not the computation time required is prohibitive. Further analysis of the well-poscdncss of the entire 
hyperbolic-elliptic system would provide a theoretical foundation for a more extensive numerical effort. 
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Appendix A: Asymptotic Expansions in Gaussian Normal Coordinates 



Quite a bit can be learned about the asymptotic behavior of the conformal factor and the regularization of the 
apparently singular terms in the evolution equations at by developing an explicit asymptotic solution. This is 
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awkward to do in our implementation of the tetrad formalism because in the elliptic equations of the Nester gauge, 
there is no local control over how the orientation of the spatial triad vectors and therefore how N a b evolves in a 
neighborhood of J^+. On the other hand, by working directly with the coordinate components of the spatial metric 
and other tensors in the CMC hypersurfaces, and by making a special choice of spatial coordinates (Gaussian normal 
coordinates based on ^+ as a two-surface embedded in the 3D conformal space [13]) in a neighborhood of J^ + , 
we can straightforwardly obtain explicit asymptotic solutions for the conformal factor and extrinsic curvature in 
terms of quantities defined from the asymptotic expansion of the conformal spatial metric as a power series in the 
conformal proper distance from J^ + . Moncrief and Rinne [l3| have done a similar asymptotic analysis in a more 
general coordinate system, but at the cost of substantial additional complexity. 

The Gaussian normal coordinates are constructed by starting from angular coordinates on the J? + 2-surface, 
which can always be chosen so that its metric is conformal to the standard unit sphere metric, with coordinates 
x A = [6, (f>) . These angular coordinates are propagated inward from ^ + along normal spacelike geodesies in the 
CMC hypersurface, and a "radial" coordinate z is defined to be equal to the proper distance from along these 
geodesies. This construction gives a conformal spatial metric of the form 

ds 2 = dz 2 + h AB dx A dx B , (Al) 

since the geodesies normal to the J? + 2-surface are perpendicular to the 2-surfaces of constant proper distance. 
Furthermore, since any metric on a 2-surface which has the topology of a two-sphere is conformal to the standard 
unit sphere metric h AB , we can decompose 

/lAB^^AB, (A2) 

defining a 2D conformal factor £ such that 

h = dot h AB = dct h AB = C A sin 2 6. (A3) 
Assuming a power series expansion of the 2D conformal metric about ^ + through order z 2 , we then have the form 



h 



AB 



Kb - 2x AB z + (x CD Xcuh A B - VVb) z 2 + O (z 3 )^ . (A4) 



Here Xab and ip AB are traceless, symmetric tensors with respect to the unit sphere metric h AB . The inverse is 



e 



+ 2x AB z + (x CD Xcuh AB + r B ) z 2 + (z 3 )] . (A5) 



It has been argued ([23], [2!|) that in a general treatment of radiation in asymptotically flat spacetimes, there should 
be polyhomogeneous terms in this expansion, i.e., zlogz and/or z 2 \ogz and higher order terms. Such terms are 
consistent with the zero-shear condition, for instance, and do not prevent the existence of .y + as a null surface in the 
conformal spacetime. As noted in Sec. HI we see no reason why such terms should be present when dealing with isolated 
systems in an astrophysical context. There may be polyhomogeneous terms present at higher orders in the expansion 
of h AB arising from the polyhomogeneous terms in the expansion of the extrinsic curvature that are associated with 
the CMC gauge choice. 

Physical interpretation of the asymptotic solutions is greatly simplified if the 3D conformal gauge is such that 
2-surfaces of constant conformal radius correspond to 2-spheres of the asymptotic flat spacetime approaching on 
a CMC hypersurface. Since the level surfaces of the 3D conformal factor approach uniformly in the conformal 
radius coordinate, this is accomplished by demanding that the intrinsic geometry of ^+ be that of a 2-sphere, i.e., 
that the 2D conformal factor £ be a constant £qi independent of angle, on J^ + , or equivalently, that the 2D scalar 
curvature of J r+ , 2 R = 2£o 2 , be uniform. Assuming that the initial data satisfy this condition, it will be preserved if 
(?t£ = 0. As a result of our asymptotic analysis, we shall see this requires a certain Dirichlct boundary condition on 
the trace K of the 3D conformal extrinsic curvature. In the context of our elliptic equations for preserving the Nester 
triad gauge, this provides the "natural" boundary condition on aK = — (3/2)Z? a £ a . The same point was noted in the 
different context of Ref. . 

The Gaussian normal coordinates are not suitable as global coordinates, since they break down somewhere in 
the interior as the normal spacelike geodesic congruence develops caustics, but they are guaranteed to be valid in 
a neighborhood of assuming a regular conformal geometry. Any breakdown of smoothness at ,f + cannot be 
attributed to the Gaussian normal choice of spatial coordinates. 
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The non-zero Christoffcl symbols for the metric of Eq. (|A1[) are just 



lh AC d z h cs = -a\ 



— ~^9 z h AB — ^ABj 



r~ A 2"pA 



where « AB are the coordinate components of the extrinsic curvature of the constant-z 2-surfaces based on the outward 
normal, which is in the minus z-direction, and where 2 r A BC are the Christoffcl symbols associated with the 2-metric 
h AB . Let k ab denote the traceless part of k ab and k the trace. Calculating the conformal 3D Ricci tensor, we get 



R z 
R z 
R A 



d z k k AB k — d z k k AB k —k , 



"|a 



(A6) 



2 6a 



d z k 



• A Z- Z- h 



! R + d z k-k 2 ) S A 



The | symbol indicates a covariant derivative in the constant-z 2-space with respect to the metric h AB . The scalar 
curvature is 



R= 2 R + 2d z k - — k - k AB k AB . 
In terms of our expansion (|A4[) of the angular metric, it is easy to verify that 

« A B = X\ + V + O (^) , 

so k A B == x A b and d z k A B = i/> A B . Also, k = (2/£)d z £. Furthermore, assuming £ = £o is uniform on J? + , 



(A7) 



2f A PA 

1 BC — 1 BC 



2 X A (bTc) - Xbc |a + (5 A (bC*c)Ko - -h BC h AD d u k 



z + Oiz 2 ) 



and 



l R = 2C 2 



( K + -K ' A | 



X"~) I z + 0(z 



(A8) 



(A9) 



The T symbol denotes a covariant derivative with respect to the unit sphere metric. 

Our definition of the conformal factor is not quite the same as in Ref . [l3| or [27j , but like them, we can solve for 
the conformal factor to third order in the expansion away from J? + . Write the expansion of f2 as 

= [1 + ci (a; A ) z + c 2 (x A ) z 2 + O (z 3 )] , 

the expansion of k as 

k = ko + k\z + O (z 2 ) , 

and substitute into the Hamiltonian constraint, Eq. (|38[) . The result is 

1 



n IT . 1. 1 

11 = Z < 1 KnZ H 

3 1 4 12 



-kl - 2ki — I K a0 K ab + k AB k 



5-AB 2 



R 



z 2 + (z 3 



(A10) 



The expansion of f2 becomes locally indeterminate, and gcnerically requires a polyhomogeneous z 3 log z term, at the 
next order. 

In the Gaussian normal coordinates, the outward normal to the constant-z 2-surfaces is just s z = — 1, s A = 0. 
Denote the traceless part of the 3D extrinsic curvature as projected into these surfaces by 



1 



(All) 



The zero-shear condition derived in Sec. I VII requires S A B = k A B = x A b- Note that x A B can be considered the 
asymptotic amplitude of the gravitational waves, and the zero-shear condition is a reflection of the fact that a 
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derivative of the spatial metric along the normal to an asymptotically null hypcrsurface becomes equivalent to a 
radial derivative of the spatial metric as the normal becomes tangent to the hypersurface. 

The momentum constraint equation |45|) allows us to determine the leading asymptotic behavior of the rest of the 
extrinsic curvature in terms of x A b- In the Gaussian normal coordinate basis, we have 



QVj [VI- 1 ki-) = Qr x d j Slk i i - 8irGn 2 ~j z 



(A12) 



First consider the z-componcnt, which through first-order in z is 



K 



a z fr 1 ^ 



1 r>A 



K 



k b a E a r - ^kK\ = ^- I 1 - ^kz 1 Q^K*. 



(A13) 



Using the zero-shear condition, the zeroth-order result is 



—K z z = —S a S b K ab = £ A B £ B A = x A bF 

30 30 BAA BA 



(A14) 



The transverse components of the momentum constraint to first-order in z arc 



—zd z { n~ l R z 
3 " 



a|b 



- kK z a = ^- I 1 - -Rz I K z A - -k\„K 
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Expand K 1 



— T B i 
a|b — ^ a|b 



\ (k z z \ , using Eqs. fl22| and (|Al4j) . to get 

V / |a 



jzd z (o-^a) + r a|b + « (d,£ B A )| B - 1~kK a - ^, B x B A = §^K Z 



(A15) 



At zeroth order, 



(A16) 



Using these zeroth-order results, we can calculate the first-order corrections. However, the equations cannot be 
satisfied by a simple power series expansion. Instead, let 



K ~ 

7^K Z Z = x A bFa + d 1 z + d' 1 z\ogz + (z 2 ) . 
Upon substitution into Eq. (|A13|) . the di's cancel, leaving d\ indeterminate locally, but 



(A17) 



d' x = -e x 



1 " AB 



(A18) 



Similarly, the coefficient of z in the expansion of (K /ZVL)K Z A is indeterminate locally, and one has an expansion of 
the form 



K ~ 

7^K Z A = x b A |b + e A z + e' A z log z + O (z 2 ) , 



(A19) 



where the coefficient of z log z is 



(A20) 



It is also possible to calculate the coefficient c' 3 of the z 3 log z term in the expansion for the conformal factor. The 
result is c' 3 = —d[/8. 

The expansions above allow us, using L'Hopital's rule, to directly evaluate at J? + the coordinate components of 
the apparently singular terms in the evolution equation for K a b, 



1 



TF 
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The result is 

S z z = \ [fofc + 2£ 2 - , 

S\ = d z (£ A B ) - V> B + ^of B - I (f D fc + 2$ - fti) <5 A b, (A21) 
<5 ,Z A = X°aTc — 2^ A ^°" 

Finally, we compute the expansion of the conformal lapse a from the CMC slicing condition, Eq. (|46l) '■ 



a = ag 



1 - \«>z + \ Qft 2 , - fti - 3x C oX D c + 2£ 2 ) z 2 + O (z 3 ) 



(A22) 



In deriving this expansion, we have assumed that So is uniform as stipulated in Sec. IVIII 

1. The Penrose Regularity Condition 

The apparently singular terms Sij are also related to the electric part of the physical Weyl tensor evaluated at J !+ . 
The Penrose regularity condition asserts that the Weyl tensor of the conformal spacetime should vanish at . The 
conformal invariance of the Weyl tensor implies that the coordinate components of the electric tensor are invariant, 

p — n \^\ v — H \^\ v — jr. ■ 

% 3 — % ^3 V ^ — % 3 ' 

A well known expression for the physical Weyl electric tensor in a vacuum spacetime (see Ref. (HI) is 

E l3 = /,', , - Ki k K kj + KKij , (A23) 

where the vacuum evolution equations have been used to eliminate time derivatives. Our assumptions about the fall- 
off of the energy-momentum tensor at imply that Eq. (|A23j) holds there. In terms of the conformal quantities, 
the traceless part of this gives 



tf k - f ~ - K ~ ~) TF r ~ ~ ~ *i TF 

En = {Rij - K t k K k] } + -Ka = - | ViVjQ + jKij | + {/.',, - Ki k K kj } 



Raising an index with the conformal metric in the Gaussian normal coordinates, applying the results in 
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Eqs. dMHSHMU), and noting the identity x A cX°b = hx C nX D cS A 



E\ = I [f D fc + 2$ ~ Hi] + I [*i - 2£ 2 - 2x C oX D c] + ^fofc = 0, (A24) 

E\ = d z (£ A B ) - V> B + \k x\ ~ I (fofc + 2£ 2 - ftr) S\ 

+ ^ A B - k x\ + \ (2£o + 2f D f c - fti) S\ - ^x C uX D c = d z (£ A B ) - h x\, (A25) 

E\ = x C aic - ^d A k + ^\a- X C a,c = 0. (A26) 

The Penrose regularity condition is not satisfied automatically once the zero-shear condition is satisfied; it requires 
the additional condition that 

9 z (S A B ) = ^oX A B . (A27) 

In the Bondi-Sachs asymptotic expansion on null hypersurfaces, the Penrose regularity condition is the vanishing 
of the traceless part of the r~ 2 term in the expansion of the metric of the two-surfaces of constant Bondi radius r 
(see [23|). In that context the condition is preserved by the Einstein equations, once imposed in the initial data, and 
therefore the same must be true for Eq. (|A27[) . 
The magnetic part of the Weyl tensor is 



Bij = Bij — VkKi^ij) 
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and using the asymptotic expansions above, one finds 

B* a =0, B Z A = 7 B* 



a z (E c B )--^x c 



with B A B = if and only if Eq. (|A27[) holds. Therefore, the magnetic part of the Weyl tensor does not give additional 
conditions. 

Given that the Weyl tensor vanishes at T + , much simpler expressions can be obtained for the apparently singular 
terms in the evolution equations than those of Eq. (|A21[) directly in terms of tetrad quantities in a general coordinate 
system, without the need to propagate the normal geodesies away from T+ to calculate K a bi etc. One just uses the 
vanishing of the tetrad components of the electric part of the Weyl tensor. The regular form of the evolution equation 

at T + for K ab is given in Eq. ((74")) . 



2. Metric Evolution Equations 

It is instructive to examine the asymptotic behavior of the evolution equations for the conformal spatial metric in 
terms of the conformal extrinsic curvature assuming the initial Gaussian normal coordinate condition is preserved. 
The general equations are 

dthj = P k d k hij + h kj d l l3 k + h lk d (3 k + 2a (k l0 + -h i:j K 



The evolution equation for h zz becomes an equation for /3 Z , 

d z (3 z = -a (k 

from which 

/3 2 =a Q 



( aK S 




It 






<SoX C dX d c + 



aK 
~3~ 



z 2 + 0(z 3 ). 



(A28) 



The condition for keeping h AZ = is 



d z p A + h AC d c /3 z 



-2aK' 



from which 



r = e 2 



h AB d B 



aK 
~6~ 



z 2 + 0(z 3 ). 



(A29) 



The boundary conditions on the shift and the lapse are from Sec. IVII1 

Decompose the evolution equations for the angular part of the metric as given by Eq. (|A4[) into trace and traceless 
parts. Defining the evolution of £ so Eq. (|A3[) is preserved, the trace in zeroth order gives 



^d t £ = y d z £ - ^aK = a Qk - ) 



(A30) 



using k = 2d z \og^. The traceless part in zeroth order, given the zero-shear condition, confirms that dth AB = 0, as 
assumed. 

The condition for keeping the intrinsic geometry of a true sphere, i.e., keeping £ uniform on T + , is d t £ = 0. 
From Eq. (|A30|) . we see this requires the boundary condition on K to be 



2aK 



(A31) 



which can be implemented in the elliptic Nestcr gauge equations of Sec. IV Al A significant consequence of this 
condition is that a — f3 z — 0(z 2 ). 
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At first order in z, and using the Gaussian- normal-coordinate shift, one finds from the trace 

D Q k = — (d t K - a Ki) = + x cd Xcd - ^-d z {aK) (A32) 

and from the traceless part 

dtX A B = a (> B - 9 2 S A B ) . (A33) 
After imposing the Penrose regularity condition, Eq. (| A27I) . 

dtx\ = 5o (i>\ - ^ox\) ■ (A34) 



2 

Next, we examine the evolution equation for the conformally rescaled extrinsic curvature, 

**' " - ^ + ^ + & + - ' *' - 2S '' + ^'1 ' 

Evaluating the right-hand side of this equation at J?+ we obtain, using Eqs. (JM! EH IATH1 IMT1 lA"22l lA"28l IA3T|) . 

8 t K z z = 0, d t K\ = 0, 8 t K\ = a Q (V' A b - 9 Z S A B ) . (A35) 

In view of the Eq. (|A33|) . this shows that the regularity conditions K z z = 0, K Z A = 0, K A B = k A B are automatically 
preserved in the time evolution. 

Now return to the question of the apparent breakdown in smoothness at J? + in the asymptotic expansion of the 

extrinsic curvature. Note that from Eq. (|A20|) . the condition for the coefficient e[ of the leading log term in K Z A to 
vanish is identical to the Penrose regularity condition, Eq. (|A27j) . However, applying the Penrose regularity condition 

to Eq. (|A18[) for the coefficient d[ of the leading log term in K z z leaves 



1/ z-2 "AB 

«1 = -?0X Tab - X 



"0 B A ~ ^K0X B A 



= -£o 2 X ab ,ab ~ ^d t (x A bX b a) , (A36) 

which does not vanish if any outgoing waves are present at J r+ . Therefore, the lack of smoothness of K z z at 
^ + cannot be prevented by a regularity condition and must be dealt with in any calculation involving emission of 
gravitational radiation. Friedrich's proofs of smoothness at J* + given smooth initial data 0, HH do not apply because 
they are based on solving a complete symmetric hyperbolic system of evolution equations. Our assumption of CMC 
hypcrsurfaccs requires an elliptic equation for the lapse, and it is this elliptic equation which generates a breakdown 
in smoothness in the presence of outgoing radiation. The elliptic equations which enforce the Nester gauge are not 
relevant because they do not affect the evolution of the coordinate metric. Our mathematical results on the presence 
of poly homogeneous terms in the asymptotic expansion of the extrinsic curvature and conformal factor, after sorting 
through differences in notation, are equivalent to those of Sec. 4 of Andersson and Chrusciel j27j . 

We do not expect the lack of smoothness at ,f + to be a serious impediment to the accuracy and/or stability 
of numerical calculations, even though in principle it could prevent rapid convergence of spectral-based numerical 
methods. 

Appendix B: Constraint Propagation 

In this appendix, we derive the constraint propagation system for our evolution scheme. This system describes the 
evolution of constraint violations which could be present in the initial data or, more generally, could be triggered by 
truncation errors in a numerical scheme. 

Since the constraint equation (|43[) is required for vanishing torsion, we should allow for a connection V with non- 
zero torsion tensor T in the constraint propagation system. An efficient method for deriving this system [4!| is based 
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on the Riemann and Bianchi identities for the curvature tensor R, 

R a ^s = (VpT a lS + T a tS T* M ) , (Bl) 

CM) (/3 7 <5) 

^ V^i? Q e7 ,5 = — R a eK sT K /3 1 , (B2) 

(/3 7 «5) (/3 7 <5) 

where X)(^ 7 5) denotes the cyclic sum over (/3j5). Expressed in terms of an arbitrary basis e a of vector fields, the 
torsion and curvature tensors are defined by 

T a 1 s = — [e 7 , e^]" — 2r a [ 7(5 ] , 
R a /3jS = e 7 (r Q ,35) — es(T a / 3 1 ) + r^r"^ — r e ^ 7 r Q e 5 — T a /3 e [e 1 , e$] e , 

respectively. If the frame e a is orthonormal, Einstein's field equations are equivalent to the requirements that the 
connection coefficients r Q( g 7 = T^p^ are antisymmetric in a/3, that the torsion tensor is zero, and that the Ricci 
tensor R a p = R 1 ' ai p is related to the stress-energy tensor T a p according to R a p = %irGT a p — AirGriaprp r^g. An 
equivalent way of stating this is the satisfaction of the equations 

T a lS = 0, A Q/ 3 75 = 0, (B3) 

where 

AafjjS = RafifS — r a ^ e T C 7 ,5 — C a pjS — S a pjS, S a pjS = 8ttG ^T]a[ 7 Ts]p ~ Vp[y T S]a ~ ^ r /a[ 7 %]/3 7 7 e '"' r eK^ ; 

and Cap^s is a tensor field satisfying the algebraic symmetries C^p^g = 0, C^p^s] = C a p 7 s = C 7 aa/3, '7 a7 Ca^ 7 <5 = 0. 
For a solution of the field equations (|B3[) . G a p 1 s is the Wcyl curvature tensor associated with the geometry. 

In order to relate Eq. (|B3|) to our evolution and constraint equations, we assume that the time-like leg eo = n 
coincides with the future-directed unit normal to the time slices E t . This implies that T 7 «5 = 0. We decompose the 
remaining components of T and A according to 

-^a 7 5 — 2_F Q [ 7 77. ( 5] -\- G a p£ 7( 5, 

A a ^ 7 a = — 4n[ Q f^][ 7 n 5 ] — e a p L V lK £ K 1 s — 2n[ a Bp] L s\s + ^a^H^ns], 

where ep 7 g = n a e a pjs an d the covariant tensor fields F, G, £, V, B and % are orthogonal to n. The left and right 
duals of A arc defined by 

(*A) Q( 9 7 5 = -e a(3 £K A eK7 5, (A*) Q( 3 7< 5 = -A a p eK e eK ~ l s, 
respectively, and induce the following transformations in the decomposition of A: 

A^*A:£^-H, V^B, B H- —V, U^-£, (B4) 
A^A* : £^B, V^H, B^-£, % n- -V. (B5) 

The trace of A is given by 

A 7 a7/3 = R aj3 + T a jS Ts 7 p - 8nG ^r a p - ^aprf 1 * 

= n a np£ 1 1 - n a ep~ l5 B~ l s + e a l5 'H 1 snp - £ a p + Vp a - (r) a p + n a np)V^ 1 . (B6) 

The Weyl tensor has a similar decomposition to that of A. However, its algebraic symmetries imply that T> = £, 
% = B, and that £ and B are symmetric and traceless. 

Our evolution system is obtained from the combinations of the field equations (|B3[) which eliminate the Weyl tensor. 
More precisely, the evolution equations (fT4l [5] [TO)) are 

= F ab B a k = -D B b k + l.o., 

= Hab = £ab ~ T^ab = ~ DoK ba + £bcdD c N c ia + l-O., 

= v a b = Bab - Hab = -D N ba - e bcd D c K da + l.o., 
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where"/. o." stands for lower order terms which do not contain derivatives of S a fc , K ab , N ab , ot6, w& or r Q( g. The 
constraint equations (fl5l [TTI [T2l nU|) are equivalent to 

= G ab B a k — —e bcd D c B d + l.o., 
= M = —T> aa = 2D a n a + l.o., 

= P a = SabcBbc = D b K a b - D a K + I.O., 

= Q a = s abc V bc = D b N ab -D a N + l.o., 

and the equation B aa — is automatically satisfied as a consequence of the hypersurface orthogonal gauge. 

Now we use the identities (|B1[ IB2[) to derive evolution equations for the constraint variables G ab , M, P a and Q a . 
For this, we first rewrite them in terms of T and A. Using the algebraic symmetries of C a p 7 s & n d the fact that 
J2(/3jS) S aM8 = 0, this yields 

(MS) (MS) 

^ V l 8(A ae7 g + C ae7 g + S ae -ys) = ^ [-r QeK A K ,9 7 5 — e fc (r ae( 9)T' v 7 (5] . (B8) 

(/3 7 5) (/3 7 5) 

We may further simplify these equations by noticing that 

^ VpS ae -y6 = -£ / 3 7 5 t e t/3 7 5 Vp>S ae y5' = £/3 7 <5 t V K (<S'*)aei.re, ^ A K /3 7 5 = £/3 7 <5 Q (A*) KCQ e , 

(/3 7< 5) (/J 7 $) 

where (S*) a| a 7 ,5 = 87rG(e 7 ,5[ Q K r / 3] K — £ a/ 3 7 ,5?7 i ''V i ,, t /3). Using the expression (|B6|) for the trace of A and the transforma- 
tion properties (|B5[) of the right dual, we find 

(A*)° £ 0e = 0, (A*) 0e be = Q b + e bcd [i cd , (A*) e ae0 = Q a , (A*) e a£h = -e abc P c - v ba + 5 a bv cc - 



As a consequence, the right-hand sides of Eqs. (|B7l IB8[) are a linear combination of the constraint fields G ab , P a , Q a 
and the evolution variables F ab , fi ab , v ab . 

Eq. (|B7|) with F a f, = yields the following identities: 

ZbcdHcd + — a a G ab = 0, (B9) 
[(Z>6 - 2n b )G ab ] D a = QaD a + G ab \D a ,D b ], (BfO) 

DoGab = —£abcPc — V ba + 8 a bV C c + K ac G cb + G ac K cb — KG a b + £acdG cb UJ d + G ac E bcd UJ d . (Bll) 

The first identity is an algebraic relation between the antisymmetric part of the evolution equation \x a b = and the 
constraint fields Q a and G ab . It can also be obtained directly by computing the antisymmetric part of Eq. Upon 
using the commutator relation s bcd D c D d = —(G ab + N ab — 8 ab N)D a , the second identity yields a differential relation 
between the constraint fields Q a and G ab . The third identity gives an evolution equation for G ab . 

In order to obtain evolution equations for the remaining constraint fields M, P a and Q a , we contract Eq. (|B8j) and 
its left dual over a/3 and then over de. Assuming that the stress energy tensor is divergence-free, this yields the twice 
contracted Bianchi identities 

V Q A 7Q 7(9 - ^A^ = *iy 5 (A*) e 5e7 - e & {T a \)T s afj + ^e s (T a ^)T 5 a7 , (B12) 

V Q (*A) 7Q 7(3 - iv fl (*A)T tf 74 = -iy J (A*) e 4e7 - e s (*r a \)T s a p + ^e s (*T a ^)T s a7 , (B13) 

where *T a p K = e a i3 1 sT lS K /2. Using the expression (|B6|) for the trace of A and its left dual, we find, in terms of the 
constraint and evolution variables, 

A 7 ° 7 o = M - Haa, A 7 % = -P b , A 7 a 7 = ~P a + Sacd^cd, A 7 al b = S a bcQc + 8 ab M - fl abl 
(*A) 7 ° 7 o = V aa , (*A) 7 % = Q b , (*A) 7 a7 o = Qa + eacdVcd, (*A) 7 a7 fc = -e abc P c + Vab- 

By virtue of Eqs. (|B9HBf Of it turns out that the zero component of Eq. (|B13[) is satisfied identically, while Eq. (|B12|) 
and the spatial components of Eq. (|B13|) give 

D M + D a P a = e bcd D b v cd + l.o., (B14) 

D P b -S bcd D c Q d + D b M = -DafXab + DbVaa+l-O. (B15) 

D Q Qb + e bcd D c P d = -D a i/ ab + D b v aa + l-o. (B16) 
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The constraint propagation system for the original evolution system, Eqs. (fT4l [9l IT0|) . is given by 
Eqs. (jBlll IB141 IB151 IB16|) where the evolution equations = v a b = are imposed. The lower order terms "l.o." 
depend linearly on the constraint fields G a b, M, Pb and Qb, but not on their derivatives. Therefore, the constraint 
propagation equations form a linear, first-order symmetric hyperbolic system, and it follows that the constraints are 
preserved by the evolution system (|141 [9l ITUf . 

In the CMC and 3D Nester gauges, the evolution equations we impose are slightly different. In order to see this, 
we decompose 

. Sab . . Sab 

Mab — fJ-ab + £abc^c + — A*, V a b — Vab + ZabcVc + 

where jl a b and v a b denote the symmetric, trace-free parts of [i a b and v a b-, respectively. Then, the equations fi a b = 
v a b = yield the evolution equations ([TBI IT?]) . Eq. (fT8|) . which yields the elliptic equation for the lapse, is equivalent to 
jj, = M. However, when writing this equation in terms of the conformally rescaled fields, the Hamiltonian constraint 
is used again to eliminate the most singular terms at J? + . This amounts to setting fj, = 0. Next, Eqs. (|T9l |2"0|) , which 
are used to derive the elliptic system for preserving the 3D Nester gauge, are obtained from v = and 2i/ a = P a , 
respectively. Using this together with 2fib = —Qb + a*a,G a b and Eq. (|B10[) . the constraint propagation system yields 

D G ab = l.o., (B17) 
D M = l.o., (B18) 

D P b -^e bcd D c Q d + D b M = l.o., (B19) 
D Qb + \e bcd D c P d = l.o., (B20) 

where Qb = Qb + o- a G a b- This system is weakly hyperbolic. However, if the Hamiltonian constraint M = is 
imposed, then the constraint evolution system is described by the Eqs. (|B17[ IB 191 IB20[) which yield a linear, first- 
order symmetric hyperbolic system with propagation speeds relative to normal observers equal to zero and one half 
the speed of light. 

Depending on the election for the shift vector in Sec. lVB| there might be additional constraints or small adjustments 
to be made in the constraint propagation system. For the spatial harmonic gauge, an additional constraint is V k = 0, 
whose propagation equation is of the form D V k = l.o. For the minimal conformal strain and distortion conditions (|61|) 
and (|B3")l , the triad vectors B a are propagated along the time evolution vector field dt instead of Do, see Eq. (|59l) . 
which corresponds to setting the combination aF a b + £bcdG ac Pd to zero instead of F a b = 0. 



Appendix C: The relation between the 3D Nester gauge condition and the 3D Dirac equation 

In this appendix, we first review the relation found in Ref. [42j between orthonormal three-frames satisfying the 3D 
Nester gauge and nowhere vanishing spinor fields satisfying the 3D Dirac equation. Then, we reconsider the equations 
in Sec. IIVCI which are responsible for preserving the 3D Nester gauge and our choice of conformal factor throughout 
evolution, and show that they can be rewritten as an inhomogeneous 3D Dirac equation. Based on these observations, 
we show how to recast these equations into second-order elliptic form. 

The Dirac equation over a three-dimensional Ricmannian manifold (£, h) with triad field B a and corresponding 
connection coefficients r a ;, c is 

ia a [Da + ^T cd a^cd) £ = mi, (CI) 
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where £ : E — > C 2 is a 5J7(2)-spinor field on E, 



1 / -i \ f 1 



1 i o J ' ° 2 = \ i o ) ' ff3 - v - 1 

are the Pauli spin matrices, Y, c d = \[a c , <Jd] and m is the mass of the field. Eq. (|C1|) is invariant with respect to local 
rotations 

£ = Ul B a = i? ab B b , (C2) 
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with U : S —> SU(2) and i? : S — > .SO^) the corresponding rotation defined by U<j a U* = a b R ba . In terms of the 
notation introduced in Sec. [Ill Eq. (|C1[) can be rewritten in the following way, 

itr (I>a - »a)£ - (m + ^n) e (C3) 

Now, suppose the frame B a satisfies the integrated 3D Nester gauge conditions ([30]) . such that n a = D a § and 
N = const. Setting m = —N/2, Eq. (|C3f) reduces to io- a Z? a (e _ *£) = with the solution £ = e*£ , £o € C 2 being a 
constant spinor. Conversely, suppose £ = (C ,*!; 1 )' 71 is a solution of Eq. (|C1|1 which is nowhere vanishing. Then, the 
rotation 

U = ^fi ~py)> p=\Z\> (C4) 

leads to the constant direction spinor solution £ = p(l, 0) of Eq. (|C3|I . It follows that the new frame B a satisfies 
n a = 5 a (logp) and N = —2m, and hence the integrated 3D Nester gauge. 

Therefore, the existence and uniqueness of a three-frame satisfying the integrated 3D Nester gauge condition (|C3[) 
with N = No = const is equivalent to the existence and uniqueness of a nowhere vanishing spinor field satisfying 
the three-dimensional Dirac equation (|C1[) with mass m = —Nq/2. The latter being a linear, elliptic equation, it 
is amenable to an analysis using standard tools from elliptic theory. Existence and uniqueness theorems for similar 
equations have been obtained in Refs. [5(| [Hi in the context of the positivity of the ADM and Bondi masses. 

While the direction of the spinor £ is related to rotations of the triad, its norm is related to conformal transformations 
thereof. Indeed, the Dirac equation (|C1[) is invariant with respect to the conformal transformation £ = Q£, B a = riB a , 
where > is a strictly positive function on S. This feature may be used in order to fix the conformal factor (up 
to a positive multiplicative constant) such that the rescaled quantities satisfy h a = 0. This corresponds to the choice 
for the conformal factor in Sec. fV\ 

Finally, consider Eqs. (ptTl I52"j) for the quantities vq = 2aD<j& and v a = aui a which are responsible for preserving 
the 3D Nester gauge and the choice of the conformal factor throughout evolution. It turns out that these equations 
arc equivalent to the inhomogeneous, 3D Dirac equation 

ia a D a i = F, F =-a(j~^A, (C5) 



for the spinor 



72 - «7i 



£=( W °- W3 ), (06) 

1 V2 — IV\ ' 



where 70 = N ab K ab + \KN + D N and j a = -e abc k bd N d c + K ab (3n b + a b ) - ^Ka a - 87rGj Q . This can be partially 
understood by considering the rotations defined in Eq. (|C4[) which map the orthonormal three-frame to one that 
satisfies the 3D Nester gauge. Parametrizing U by its rotation angle Q and axis e a , e a e a = 1, we have 



U = exp ( Yi ea ° a ) = cos ( "2 J ~ lsm \ ~2 ' ea<7a ' 



Hence, the corresponding spinor field is 



I cos(f) -isin(|)e 3 
(e 2 -i ei ) S in(|) 



For an infinitesimal rotation and change p = po + Sp in the conformal factor, like the one required to preserve the 3D 
Nester gauge and our condition on the conformal factor in an infinitesimal time step, this yields precisely a spinor of 
the form (|C6[) with vq = Sp and v a = poSOe a /2. 

Eq. (|C5j) constitutes an elliptic system of equations for the fields vq and v b which can be solved subject to appropriate 
boundary conditions. A second-order elliptic system can be obtained from this by the ansatz £ = ia b D b r\, with rj a 
potential spinor field. Introduced into Eq. (|C5|) this yields 

- DaDai] + iN ab a b D a r] - iNa a D a j 1 = F, (C7) 

N bc D b — ND C . In terms of the potentials ip and £ a 



introduced in Sec. IV Al the potential spinor is given by 



C2 - *Ci 



(C8) 
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